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SUMMARY 

The flexible parallel link mechanism is designed for increased rigidity to sus- 
tain the buckling when it carries a heavy payload. Compared to a one link flexible 
manipulator, a two link flexible manipulator, especially the flexible parallel mecha- 
nism. has more complicated characteristics in dynamics and control. The objective 
of this research is the theoretical analysis and the experimental verification of dy- 
namics and control of a two link flexible manipulator with a flexible parallel link 
mechanism. 

Nonlinear equations of motion of the lightweight manipulator are derived by 
the Lagrangian method in symbolic form to better understand the structure of the 
dynamic model. The resulting equations of motion have a structure which is useful 
to reduce the number of terms calculated, to check correctness, or to extend the 
model to higher order. A manipulator with a flexible parallel link mechanism is a 
constrained dynamic system whose equations are sensitive to numerical integration 
error. This constrained system is solved using singular value decomposition of the 
constraint Jacobian matrix. Singular value decomposition is a stable algorithm for 
the dynamic analysis of a constrained system. Elastic motion is expressed by the 
assumed mode method. Mode shape functions of each link are chosen using the load 
interfaced component mode synthesis. The discrepancies between the analytical 
model and the experiment are explained using a simplified and a detailed finite 
element model. The step response of the analytical model and the TREETOPS 
model match each other well. The nonlinear dynamics is studied using a sinusoidal 
excitation. The nonlinear dynamics due to the flexibility is significant. However, 
the nonlinearity of RALF (Robotic Arm, Large and Flexible) is not fully studied 
experimentally due to the speed limitation of the hydraulic cylinder. The actuator 



xvii 

dynamic effect on a flexible robot was investigated. The effects are explained by 
the root loci and the Bode plot theoretically and experimentally. For the base 
performance for the advanced control scheme, a simple decoupled feedback scheme 
is applied. 
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CHAPTER I 


Introduction 


1.1 Motivation 

A large two degree of freedom flexible manipulator designated RALF (Robotic 
Arm, Large and Flexible) as shown Fig. 1.1 has been constructed in the Flexible 
Automation Laboratory at Georgia Tech. The structure consists of two ten foot long 
links made of aluminum tubing actuated by hydraulic cylinders. The upper link is 
driven by a parallel link mechanism. This large manipulator cam reach farther than 
a traditional robot. Such a configuration would be useful for material handling, 
for welding, or for ultrasonic inspection of a large structure such as an airframe . 
Using a lightweight material, it is possible to construct a large manipulator with 
low power consumption and high speed operation. However, the control law must 
take into account the reduction of the structural vibrations due to the distributed 
flexibilities of links. 

In a conventional serial link mechanism, the upper link is driven by a motor 
attached at the tip of the lower link. The weight of the second motor must be 
carried by the first motor which increases the torque and power the first motor 
must provide. Especially in lightweight arms, the weight and reaction torque of the 
second motor increases the structural vibration of the lower link. The structural 
vibration makes the lower link fatigue more quickly. To reduce these effects, the 
rigidity of the lower link should be increased. One possible remedy is to include 
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a parallel drive mechanism. The advantages of parallel drive mechanisms are high 
rigidity, high loading capacities, and lower interaction between joints. 

Although a variety of research on serial link mechanisms can be found in the 
robotics literature, the studies of closed-loop chains such as the parallel mechanism 
are rare and mostly concerned with rigid manipulators. In a parallel mechanism 
with rigid links, there is a simple relationship between generalized coordinates. 
However, in a parallel mechanism with flexible links, this relationship is expressed 
in the form of a complicated nonlinear algebraic constraint equation. To solve the 
dynamics of a flexible parallel link mechanism, a mixed set of differential equations 
of motion and nonlinear algebraic equations must be solved simultaneously. Because 
the numerical integration of these equations is sensitive to numerical error, a stable 
integration algorithm is required. 

One link flexible manipulators have been studied in both theoretical and ex- 
perimental aspects recently. However, a different approach is needed in dynamic 
analysis and control of two link flexible manipulators. 

First, two link flexible manipulators have highly nonlinear interaction terms be- 
tween links. Furthermore, dynamic equations of motion of flexible manipulators are 
more complicated than those of rigid manipulators. The number of equations of 
motion increases a s the number of modes to be included increases. It is difficult to 
understand the effect of flexible motion on rigid motion via recursive forms of the 
equations of motion for a multi-link arm even if it were efficient to derive a inverse 
dynamics. On the other hand, a closed form of the equations of motion is useful in 
understanding the characteristics of model parameters. However, the equations re- 
sulting from existing closed forms Eire too complex to serve this purpose. Therefore, 
a method which is structurally well organized and computationally efficient must 


be developed. 
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Secondly, mode shape functions of each link are not simple to derive due to 
coupling between links. The common approach to determine modes has been the 
finite element method. Even though the finite element method is a systematic 
modeling technique for complex mechanisms, computational efficiency is lost due 
to the fact that a large number of elements are needed to obtain appropriate and 
accurate results. Because of its use in control design, the resulting model should be 
simple enough to render the analysis at hand tractable while retaining the significant 
features of the original structure. For simple structural shapes, a continuous system 
approach using assumed modes seems to be an easier way. Using appropriate modes, 
a lower order model can be obtained. 

Third, the control of a flexible manipulator has different characteristics com- 
pared to the control of a rigid manipulator. Even though control inputs are applied 
only at the joint, rigid body motion and elastic body motion have to be controlled 
simultaneously. The control law constructed on the basis of a reduced order model 
must have robustness for truncation error and for uncertainty of its parameters. 
The control bandwidth with joint feedback alone is below the lowest vibration fre- 
quency. Low control bandwidth means a long settling time and poor disturbance 
rejection to external forces. Increased control bandwidth can be obtained by feed- 
back of additional states such as strain and tip position which provide information 
on the elastic motion of the manipulator. The approximate mode shapes are neces- 
sary to get these additional states. Whereas the one link arm is essentially linear, 
the two link arm is strongly nonlinear and its dynamics change according to its 
configuration. 
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1.2 Review of Previous Work 

This review consists of several topics which are mentioned in the previous section 
on motivation. 

1.2.1 Derivation of Equations of Motion 

One of the primary concerns in manipulator dynamics is computational effi- 
ciency. For the efficient form of the manipulator dynamic equations, various re- 
cursive formulations for rigid manipulators using Lagrangian [26], Newton - Euler 
[42], or Kane’s method [34], have been proposed. For flexible manipulators, Book 
used the method of homogeneous transformation matrices. He first considered small 
linear motions of a massless elastic chain [7] and later considered distributed mass 
and elasticity [10]. When the recursive formulation is used, the structure of the dy- 
namic model which is quite useful in providing insight for designing the controller 
is destroyed. To overcome this problem, several programs for rigid manipulators 
have been developed to derive the equations of motion in symbolic form [13,52]. 
Symbolic formulation has the advantage of allowing the identification of the dis- 
tinct components of the model. For flexible manipulators, Maizza-Neto [45] derived 
symbolically the equations of motion of a two link flexible manipulator by hand. 
A systematic method to symbolically derive the nonlinear dynamic equations of 
multi- link flexible manipulators was presented by Cetinkunt [14]. However, he did 
not explore the structure of the terms in the flexible manipulator model. The con- 
ceptual framework leads to guidelines for simplifying robot dynamics. The physical 
interpretations and structural characteristics of the Lagrangian dynamic model for 
rigid robots were derived by Tourassis and Neuman [71,72]. The inertia matrix is 
induced by the masses and center of mass of links. In turn, the centrifugal and Cori- 
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olis coefficients are derived from an inertia matrix through the Christoffel symbol. 
But, the method of deriving mass matrices is not efficient for a flexible manipula- 
tor. Asada [2] presented a method which uses the Jacobian matrix to derive the 
mass and gravity matrices. His method is found in this thesis to be ver\ efficient in 
modelling a flexible manipulator. Low [41] used the Jacobian matrix in deriving the 
equations of motion of a flexible manipulator. However, the structure of centrifugal 
and Coriolis force is still complicated and hard to understand. 

1.2.2 Mode Shape Function 

A flexible manipulator that undergoes large rigid body rotation can be modelled 
in terms of either an unconstrained modes method where the entire body vibrates, or 
a constrained modes method where the joint is held motionless. The unconstrained 
modes approach is more rigorous. However, the calculation of the natural frequen- 
cies and the mode shapes is relatively complicated and it is difficult to extend this 
approach to multi-link arms. On the other hand, the constrained modes method is 
an approximate solution of the unconstrained modes method when the beam to hub 
inertia is small enough. Because with the constrained modes method it is simple to 
derive the natural frequencies and the mode shapes, it can be applied to multi-link 
arms. The unconstrained modes and constrained modes methods were compared 
by Barbieri [3]. Hughes [29] and Hablani [23] compared the accuracy of these two 
methods by modal identities and completeness indexes respectively. Schmitz [62] 
and Hastings [25] experimentally verified their model with unconstrained modes 
method and constrained modes method respectively. In the case of multi-link arms, 
the boundary conditions at the joint are not clear because there is coupling between 
links. In most studies of multi-link flexible arms, the boundary conditions of each 
link are assumed to be clamped -free. The accuracy of this assumption has not been 
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yet experimentally confirmed. 

The reduced order model for design of controllers can be obtained by truncating 
the modes with frequencies higher than the actuator bandwidth. This is reasonable 
because high frequency modes generally have small amplitudes and can be regarded 
els disturbances. Furthermore, actuators which are not excited by the high frequency 
modes act as lowpass filters . Hughes [30] and Skelton [68] investigated the criteria 
of model order reduction. Hastings [25] and Schmitz [62] got good experimental 
results using only two or three modes. Tsujisawa [74] derived a reduced order 
model for RALF by the modal cost analysis method. He suggested that using the 
first two modes of each link is optimal from the control point of view. 

The elastic deformation of links can be expressed by a set of admissible functions. 
The selection of admissible function sets is not unique. The finite element method 
(FEM) uses simple admissible functions. FEM can easily deal with complicated 
structures with complex boundary conditions. However, a large number of elements 
axe required to obtain accurate modes and frequencies. To reduce the dimension of 
the model, a component mode synthesis (or substructure modal synthesis) has been 
developed. A complete structure can be regarded as an assembly of substructures. 
The lower frequencies and corresponding mode shape functions of the complete 
structure can be obtained by synthesizing a truncated set of modal properties of 
the substructure. Three basic methods have been developed depending on whether 
its interface is held fixed [19], free [27], or loaded [6]. Sunada [70] and Shabana [66] 
obtained a nonlinear finite element model of a flexible structure for a large rigid 
body motion by using the component mode synthesis. However, these programs 
still execute too slow for a real time control. 
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1.2.3 Constrained Dynamic System Analysis 

The dynamics of robot manipulators containing general closed chains in their 
structures has been studied very little. The parallelogram, a special case of the 
closed chain, has been well studied by Asada [1]. Chung [18] derived the equations 
of motion of RALF, but he assumed that the actuating link and the lower link 
had the same mode shapes so that the constraints did not need to be considered. 
Megahed [48] and Luh [43] derived the equations of motion of rigid robots with 
closed chains by Lagrangian and Newton-Euler methods respectively. However, 
a computationally stable algorithm is required because the model is sensitive to 
numerical error. 

There are two conceptual approaches to solve the constrained equations of mo- 
tion. One approach solves the equations of motion simultaneously with the con- 
straint equations. The other approach uses a reduction method that eliminates the 
constraint forces explicitly from the equations of motion. It is difficult to solve a 
mixed set of differential equations of motion (dimension n) and nonlinear algebraic 
equations of kinematic constraint (dimension m). For this purpose, Nikravesh [54] 
reviewed three integration algorithms. These are the direct integration method, the 
constraint violation stabilization method and the generalized coordinate partition- 
ing method. The direct integration method [48,43] converts the algebraic equations 
to second order differential equations, then solves these equations with equations 
of motion. This method results in constraint violations because, as time progress, 
the integration numerical error accumulates. The constraint violation stabilization 
method [4] introduces constraint violations as a feedback term to correct the vio- 
lations in the next integration step. One difficulty of this method is the selection 
of proper feedback gains. These two methods are sensitive to initial conditions on 
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the system. In the third method [76], the generalized coordinates are partitioned 
into independent (dimension n - m ) and dependent (dimension m) sets. Numer- 
ical integration is solves for the independent generalized coordinates. The choice 
of the correct initial condition is not critical and the dimension of the equations of 
motion is minimum. When generalized coordinates are partitioned, an important 
consideration is the choice of independent generalized coordinates. An arbitrary 
selection of independent generalized coordinates often results in ill - conditioned 
matrices. Wehage [76] identified independent generalized coordinates by using LL 
partitioning of the constraint Jacobian matrix. Maui [46] and Singh [67] used sin- 
gular value decomposition, and Kim [37] used QR decomposition. A unique and 
accurate independent set of generalized coordinates can be obtained from the last 
two methods. Singular value decomposition is a more robust algorithm than the 
QR method and has a variety of applications in linear control systems. 

1*2,4 Motion Control of a Flexible Robot 

A colocated control scheme, in which sensors are colocated with actuators, has 
been used for most flexible manipulator control. With a PD control for each joint 
(Independent Joint Control), Book [7] showed that the maximum closed loop band- 
width of rigid body motion is a half of the system first natural frequency. If inter- 
joint feedback terms are included between actuators (Generalized Rigid Control), 
the closed loop bandwidth can be increased up to the system first natural frequency. 
Golla et al. [21] showed that IJC also provided control as good as GRC. 

A noncolocated control, where the tip position was measured by an optical sen- 
sor, was implemented by Schmitz [62]. Because tip position feedback alone creates 
a nonminimum phase system, the achievable gain is limited. However, improved re- 
sponse, two times faster than the first natural cantilever frequency, was obtained by 
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the combined feedback of tip position, link strain, and hub rate. Better disturbance 
rejection to end position external forces than with joint feedback alone was also 
obtained . However, the response time is limited by the inherent wave propagation 
delay for the beam. A precise dynamic model and a sophisticated control law are 
required for satisfactory response. 

Control of flexible robots is characterized by control concepts similar to those 
being used in the control of rigid body robots. Four major methods have been devel- 
oped: linear feedback control, computed torque control, adaptive control, singular 
perturbation control. 

For a one link flexible manipulator, Hastings [25] implemented a Linear Quadratic 
Regulator with full state feedback. The time varying modal amplitudes were re- 
constructed from strain measurements. A reduced order observer was utilized to 
obtain estimates of the modal velocities from the reconstructed modal amplitudes. 
Sakawa [61] also used LQR. Schmitz [62] used a Linear Quadratic Gaussian regu- 
lator and its reduced order compensator. These compensators were experimentally 
verified regarding robustness and disturbance rejection. Krishnan and Vidyasagar 
[40] experimentally showed that the performance obtained using the bounded input 
H 2 optimal controller was better than that obtained using a discrete time LQG con- 
troller. For a two link flexible manipulator, Maizza-Neto (45] discussed the use of 
a pole placement algorithm to obtain full state feedback gains. Full state feedback 
showed high sensitivities to parameter perturbations and higher torque require- 
ments. Ower and Van De Vegte [56] applied a classical design technique using Bode 
plots to a transfer function model of a two link flexible manipulator. The last two 
methods were not experimentally verified. 

Unlike one link flexible manipulators, the two link flexible manipulator is strongly 
nonlinear due to interactions between links. Bavo [5] used a computed torque 
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method for trajectory control. The torque at each joint can be found in the fre- 
quency domain by means of an iteration procedure. However, the iteration aspect 
prevents this method from being used in real time. Pfeiffer et al. [58] developed a 
multistage control scheme. First, the rigid body motion is controlled by a conven- 
tional computed torque method. Second, the elastic deviations from the reference 
path are quasi-statically corrected by modifying the reference path. Third, the 
remaining elastic vibrations are actively damped by strain feedback for each link. 
Schutter et al. [63] presented a nonlinear feedback which linearized the rigid body 
dynamics, followed by linear feedback of the full state. The last two methods, like 
other computed torque methods, lead to a computationally intensive controller, 
sensitive to model parameters. 

Various adaptive concepts are categorized into two types; MRAC(Model Ref- 
erence Adaptive Control) and self-tuning control. Siciliano, Yuan and Book [65] 
proposed a full state type MRAC for one link flexible arm. Meldrum and Balas [47] 
used direct adaptive type MRAC, in which the controller parameters are adjusted 
with only the plant output and input signals. Yuh [83] and Yang and Gibson [80] 
presented an indirect adaptive control approach based on an identified linear predic- 
tion model of the plant. Rovner [60] developed an adaptive algorithm based on the 
self- tuning regulator concept for the noncolocated control case and experimentally 
proved its performance. Nelson and Mitra [53] and Yurkovich and Pacheco [84] pre- 
sented load estimation and load adaptive control. Most researches and experiments 
of adaptive control are performed for one link flexible arms because an assumption 
that during the adaptation process the elements of the linearized system remains 
constant is well satisfied for one link flexible arms. However, in multi-link flexible 
arms, the controller has to be applicable to rapidly time varying nonlinear cou- 
plings. Cetinkunt [15] applied Adaptive Model Following Control (AMFC) based 
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on the generalized inertia matrix for a two link flexible robot. This controller relaxes 
some of the restrictive assumptions made by previous AMFC design procedures so 
that the use of the AMFC techniques in high speed manipulators becomes possible. 
Yuan [82] applied a robust controller based on MRAC for RALF. 

Under the assumption that only small elastic deviation from rigid body motion 
occurs, the decoupled control of rigid and elastic coordinates can be considered. 
Such a concept leads to a two stage control scheme, a slow control for the rigid 
body motion and a fast control for the elastic motion. Siciliano and Book [64] 
applied a singular perturbation method to a one link flexible arm. 

1.3 Thesis Outline 

In chapter 2, nonlinear equations of motion of RALF are derived by the La- 
grangian method in symbolic form. The derivation procedures for the equations of 
motion are described in detail. 

In chapter 3, the nonlinear algebraic constraint equations of RALF are derived. 
A mixed set of differential equations and algebraic equations are solved using sin- 
gular value decomposition. The main concept of this method is explained using a 
simple example. 

In chapter 4, the method for determining mode shape functions of each link 
is described. The loaded interface component mode synthesis is used to find the 
proper mode shape functions. The discrepancies between the predicted frequencies 
and the measured frequencies axe explained by the finite element models. 

In chapter 5, first, the time response of the analytical model is compared to that 
of the TREETOPS model for the model verification. Secondly, the nonlinear dy- 
namics of RALF axe studied by sinusoidal excitation. Third, the effect of hydraulic 
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actuator dynamics on flexible arm dynamics is discussed. 

In chapter 6, a decentralized control scheme using cylinder position and strain 
of each link is applied to check the characteristics of the control of a two link flexible 
robot. 

1.4 Contributions 

The major contribution of this thesis is the theoretical analysis and the experi- 
mental verification of dynamics and control of two link flexible robots with flexible 
parallel link mechanisms. 

Detailed descriptions of contributions axe as follows: 

First, nonlinear equations of motion of RALF are derived in symbolic form 
systematically and efficiently. The resulting equations of motion have a structure 
which is useful to reduce the number of terms calculated, to check correctness, or 
to extend the model to higher order. 

Secondly, the dynamics of a closed kinematic chain system with a flexible parallel 
link mechanism is solved without any significant constraint violation. 

Third, the proper mode shape functions of each link of RALF are determined 
using component mode synthesis. It is verified that component mode synthesis 
provides rigorous boundary conditions for modal data of components. 

Fouth, nonlinear dynamics of RALF is verified using a sinusoidal excitation. 
The degree of nonlinearity of RALF is observed from the power spectra of the tip 
acceleration. 

Fifth, the effect of actuator dynamics on the flexible robot dynamics is studied. 
The hydraulic actuator dynamics is modeled and compared with the dynamics of 
the electric D.C. motor. The difference in velocity feedback for the two actuators 
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is shown to result in significantly different joint behavior for a flexible structure. 

Sixth, a decentralized control algorithm using cylinder position and beam strain 
has been demonstrated for the control of RALF. It is shown through experiments 
that position feedback through a lag compensator and strain feedback through an 
nonminimum allpass filter yield a good trajectory following and beam vibration 


suppression respectively. 
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CHAPTER II 

Derivation of Equations of Motion 

2.1 Introduction 

In this chapter, a Lagrangian method is used to derive the equations of motion 
for a flexible manipulator with parallel link mechanisms. Because flexible manipu- 
lator dynamics is more complicated than rigid manipulator dynamics, the amount 
of computation increases and there is a possibility of ending up with incorrect 
equations. Therefore, an efficient and systematic derivation method is required to 
reduce these problems. Furthermore, it is desirable to simplify elements of the mass 
matrices and the centrifugal and Coriolis forces for a real time control. 

2.2 Description of Structure 

The structure as shown in Fig. 2.1 consists of lower,, upper, connecting, and 
actuating links. Each link is connected to another by a pin. The upper link is driven 
by a parallel link mechanism. Motion is restricted to the vertical plane. The joint 
parts of the lower and the upper link axe stiffened by increasing the cross sectional 
area. Detailed structural data are given in Wilson [79]. In deriving equations of 
motion using the Lagrangian method, the coordinates involved in a closed kinematic 
chain are not independent. In the flexible parallel link mechanism, the relationship 
among coordinates axe described by nonlinear algebraic equations. Therefore, a 
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constraint dynamic analysis is required to solve nonlinear algebraic equations and 
differential equations of motion simultaneously. To derive the equations of motion 
of this closed kinematic chain system, one joint of the parallel link mechanism is 
virtually cut to form an open tree structure as shown in Fig 2.2. The unknown 
constraint force is applied at the virtually cut joint. 

2.3 Equations of Motion 

The open tree structure is regarded as an assembly of two serial link manipu- 
lators - the lower and upper link part plus the connecting and actuating link part. 
The equations of motion of each part can be written symbolically. 

E + E + E E c,k&mk + a. g = r, ( 2 . 1 ) 

j=\ j = l j=i fc=i 

where q 3 is an element of the generalized coordinate vector, M t] is an element of 
the generalized mass matrix, A' tJ is an element of the elastic stiffness matrix, CAi) 
is an element of the velocity coupling matrix which is the coefficient matrix of 
Coriolis and centrifugal force, G, is an element of the gravity force matrix, g is the 
gravitational acceleration vector, r, is an element of the generalized force vector. 

The mass matrix and the gravity force matrix can be derived using the Jacobian 
matrix as shown in Appendix A. 

Mtj-Zf’ JJJ t p r A r ix r (2.2) 

P= 1 J0 

where b is the number of links. 7 P , p p , A p , and l v are the Jacobian matrix, the 
density, the area, and the length of link p respectively. 

g, = E [’ J';'p r A r dx p 


(2.3) 
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where J ^ is the i th column of Jacobian matrix J p . The velocity coupling matrix 


can be derived from the mass matrix. 


Cjkd) = =( 


1 .dMn , dM ik dM jk 


+ 


) 


2 V dq k ' dq : dq, 

The stiffness matrix is related to the mode shape function ti' :j . 

d\- t 


Ui = f' E,I,( 

JQ 


'U \2 




)‘dx t 


(2.4) 


(2.5) 


Next, equations of motion of two parts are combined using unknown constraint 
force for the complete equations of motion. Equations of motion of the closed 
kinematic chain system can be written svmbolicallv as follows. 


H K U9j + H C )k(i)Wk 

i - 1 j=i j=i fc=i 


(2.6) 


+ Gig + — t, 

k=\ 

where m is the number of the constraint equations, A* is element of the unknown 
constraint force vector, ($ ? )jti is element of the constraint Jacobian matrix which is 
derived by differentiating the constraint equations (2.7) with respect to time. 

(2.7) 
(2.S) 


*(9) = 0 

d*(q) d*(q) dq 


dt dq dt 0 

In the following sections 2.4 and 2.5, the derivation procedure of the equations 

of motion of each part is described in detail. The constrained dynamic analysis will 

be described in Chapter III. 


2.4 Lower and Upper Link 


2.4.1 Mass Matrices and Gravity Force Vectors 


Deformed position vectors of each link in Fig. 2.3a and 2.3b are described as 


follows: 
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fi = (x\cos@\ — u-isin&i )i + ( Xisin&i + UjCos9i)j (2.9) 

f*2 = [liCosOi — ui e 5in^ + X2Cos(9\ + $ 2 ) — uisin(9 1 + 0j)]i 

4- [/ism#i + u\ e cos9i 4- X2sin(9\ + ^ 2 ) + U2Cos(9i 4- ^2 )]j (2.10) 

where i and j Eire unit vectors Eilong the inertial frame, Xo and Vo- The elastic 
deformation, u*, can be expressed by finite series of mode shape functions which 
satisfy assumed boundary conditions multiplied by time dependent modal coordi- 
nates. Suppose that the amplitude of the higher modes is relatively small compared 
with the first mode, two modes per link Eire considered in this model. 

ui(xi,t) = ^n(*i)6i (<) + 012(*l)£l2(<) (2-11) 

U 2 ( x 2 , t ) = ^2l(*2)6l( < ) + V»22(^2)62(<) (2-12) 

The elastic displacement of the end point is 


U U = Ui(/i,t) 


(2.13) 


Velocity vectors cire related to generalized coordinates by the Jacobian matrices [2]: 


r*i = Jiqu 


(2.14) 


f 2 = J2912 


(2.15) 


where the generalized coordinates are 


512 = {0l>02,fu>£l2>f21>f22}^ 


(2.16) 


and the Jacobian matrices Eire 


Ji = 


—u\C\ — ii S\ 0 — i/^Ci 9 9 

— xi\S\ -f- x\Ci 9 rinCi 1P12C1 9 9 


(2.17) 
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J 2 = 


— l\S\ — U\ e C\ — U;Cn ~ X2S12 — U 2 C X2 — X 2 S \ 2 
+ l\Ci — Uj e Si — U 2 -S'i 2 + X2C12 — U 2 * 5 i 2 “H X 2 C \ 2 


~ty\uS\ ~ ^\2eS\ ~ 021-S'i2 — U> 22 S\ 2 

iblleCl ^\2eC\ 021^12 ^ 22 C\ 2 

where co.s(<9 t + 9j) and sin(9 t + 9,) are expressed as C x: and Sij for convenience. The 
Jacobian matrices, J\ and J 2 , can be easily derived from the position vector using 
the MJac function of SMP (Symbolic Manipulation Program) [69]. Using these 
Jacobian matrices, mass matrices and gravity force matrices are calculated by the 
following equations: 


(2.13) 


Mi, = / J x J\P\A x dx x + f J 2 r J 2 p 2 A 2 dx 2 (2.19) 

Jo Jo 

G x = [ h J x [2,l}p x A l dx 1 + f 2 J 2 \2,l}p 2 A 2 dx 2 
Jo Jo 

G 2 = f 7 J 2 [2 1 2]p 2 A 2 dx 2 (2.20) 

Jo 

J p [2, i] is the second row and the ith col umn of J p selected since the gravity is acting 
in the negative direction of Yq. The gravitational potential energy change due to 
the link deformation is assumed to be negligible. 

Elements of mass matrices and gravity forces are: 


ri\ 

-Mn = / (*i + u\)p\A\dxi 

Jo 

+ / V? + + «2 + *2 

Jo 

+ 2 (/!X 2 C 2 — S*i + U le U 2 C 2 + ^\e x 2^2))P7^2 ( i x 2 

J r^i 

i ( x\ + u\ -I- h^C 2 — l\U 2 S 2 + U\ e U 2 C 2 + Ux e X 2 S 2 )p 2 A 2 dx 2 
0 

ft i rh 

Mi 3 = / XitPuP\A\dxi + t/’ne / (k + x 2 C 2 - u 2 S 2 )p 2 A 2 dx 2 
Jo Jo 
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Mi 


14 


= f X\4>\2p\A\dX\ + 0i2e f {h + x 2^2 V>2^2)P2^ 

JO 

rh 

Mis = / ^21 ( x 2 + h^2 + U\ c S2)P2^2d^2 

Jo 

rh 

M\ 6 = / 022(^2 + ^ 1^2 + Ui c 52 )P 2 - 4 . 2^2 

Jo 


2 ^X 2 


A 2 o o 

A/22 = / 0 2 + ^2)^2 ^2 dX 2 

Jo 

rh , 

M 2 3 = V^lle / ( £ 2 ^2 - U 2 5 2 )P 2 ^ 2^2 

Jo 

Wj - 

M 24 = V> 12 e / (Z2C2 “ U 2 52 )P 2 ^ 2 «X 2 

Jo 

rh 

M 25 = / xi^-npiAidxi 

Jo 

Mis = / x 2 ^ 22 p 2 A 2 dx 2 

Jo 

M33 = / ^nPiAidxi + rplu f piAidxi 
Jo Jo 

rh f^ 2 

M34 = f 0 H^ 12 Pl-Al^ x l + ^Ue 012 e / P2-^2^ x 2 

Jo ' /0 

rh 

M35 = rpueCi / V> 2 1 piAidxi 

Jo 

rh 

M36 — V'ne^ / rpuPiAidxi 
Jo 


l W] 

M 44 = / 1 ^uP'Aidxi +1>lu I PiA-idxi 
Jo J ° 

M 45 = TpUeCl f faiplAidli 
Jo 

M 46 = *I>12 eC-i f ifrllPlAldXi 

Jo 


fh 

Mss = / 'p\\P'iAidx2 

Jo 


( 2 . 21 ) 
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r h 

M$6 = / ^ 21 022 /> 2 - 4 . 2^2 

JO 


Mi 


66 


f 7 2 

= / 022 ^ 2 - 
JO 


A 2 dx 2 


Gi = [ (xiCi - UaSiVi^! Jxi + f {l\Ci - U\S\ + X 2 C\2 - 'U’2S\2)plMdx2 
Jo Jo 

G 2 = / (^2^12 “ W2*Si2)p2-42<iX2 (2.22) 

Jo 


where 




SMP files of equations (2.9 - 2.22) are included in Appendix D. 
The integral in the above elements axe defined as follows. 


TTX{ — / 

Jo 

(2.23) 

fix 

miU c = / XiP\A\dxi 
Jo 

(2.24) 

Wi 

Di = / x'piAidxi 
Jo 

(2.25) 

rit 

LMij = / ^\j(^Xi)pxAxdxi 

Jo 

(2.26) 

AMij = / Xi\l>ij(xi)piAidxi 

Jo 

(2.27) 

NM tJ = [' tfjix^piAidxi 

(2.28) 


where l ie is center of mass of link i. The first three terms are parameters which are 
related to a rigid motion. These are called zeroth, first, and second moments of 
inertia respectively. On the other hand, the last three terms are parameters which 
are related to a flexible motion. LMij and AMij are called the modal momentum 
coefficients and the modal angular momentum coefficients respectively [29,30]. The 
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physical meaning of these terms is not easy to explain, 
following properties [29.30] . 

£ LM l = m ' 

j - 1 


However, these have the 
(2.29) 


£ LMijAMij = mJi C 

j=i 

;=i 

NMij is used for the normalization of mode shape functions. 

Using these six inertia parameters and the orthogonality of modes and linearizing 
about zero deflection, the linearized mass matrices and gravity vectors are rewritten 
as follows: 


(2.30) 

(2.31) 


Mn = Di + Di 4- + 2 / 2 ^ 2 ) 

■Ml 2 = D 2 + 77*2^1 ^20^2 
M\z = AMn + tywcTniih + hcCi) 
Mu = AM 12 + ^\2e m 2(h + hcCz) 
Mi 5 = AM 21 + LM 21 hCi 
M16 = AM22 + LM22 I1C2 

M22 = &2 

M23 = l/’lle m 2^2cC2 
M 24 = 1pUe^2hcC2 
M 25 = AA /21 


M 26 = AM 22 


(2.32) 
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Mzz = N Mu + m 2^\\ c 

M 3 4 = 0 

A /35 = ^lic -^A/ 21 C 2 

A/36 — 0ncIA'/„C 2 

A/u = jViV /12 + m 2 0i 2e 

A/45 = 01 2c ^ A/21 ^2 

A/46 = 012c A/22 C*2 

M 55 = A/ A /21 
A /56 = 0 


A/s6 — A" A /22 


Gi = mi/i c Ci + ^ 2 (^ 1 ^ + ^ 20 ^ 12 ) 

G 2 = n22/2cC > i2 (2.33) 

Finally, the elastic stiffness matrix, ivYp can be derived by the partial differenti- 
ation of the elastic potential energy, V : 

v = 5 f Ei ^ d *' + 5 C **%%?*. (2 - 34) 

K„ = £ ElS^Vdzi ( i,j = 1,2) (2.35) 

where E is Young’s modulus of elasticity, and Ii is the area moment of inertia of 


link i. 
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2.4.2 Payload 

When masses are loaded at the end of the lower link and the upper link, position 
vectors of the payload are derived from (2.9) and (2.10) by substituting x, into /,* 
and Ui into u ie . 



f* le = (l\cos8\ — U\ e sin6i)i 4- 4- u\ e cos8i)j 

(2.36) 


f* 2 e = [Zicoa^i — Ui e sin9i 4- l 2 cos(8i + 0?) — u 2e sin($ 1 4- ^ 2 )] 1 



4- 4- u\ e cos9\ 4- l 2 sin(9i 4- # 2 ) 4- u 2e cos{8\ 4- # 2 )!/ 

(2.37) 

where 

Ui e = 

(2.38) 


U2e = U2 (l 2 ,t) 

(2.39) 

Velocity vectors are expressed by Jacobian matrices: 



r le — J le?12 

(2.40) 


>~2e = «^2e<?12 

(2.41) 


where 

-liieC! - /iS! 0 -0n.5i -iha.Ci 0 0 

Ju = 

— 0 1pl2eC\ 0 0 

— /j5i — UieCi — ti2e^ r i2 ~ h$12 ~ u 2eCi2 ~ ^2*?i2 

J 2 e — 

+/lCi — Ui c 5i — ^26^12 + ^2^*12 “^2c‘5'l2 + ^2^12 

“ ^21c^l2 ~“022e*Sl2 

ll>U eC\ lt>\2eC\ *l>2uCi2 *I>22 cC\2 

Mass matrices and gravity force vectors due to the payload, m lp and m 2p , are 

expressed by Jacobian matrices. 

Mfj = rri\ p jJ t Ju + m 2p^2^2e 


(2.43) 



(2.44) 
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G\ = m lp J le { 2, 1] 4- m 2 pj 2e [2, 1] (2.45) 

G p 2 = m 2p J 2e [ 2,2] 

Following the procedure of the previous section, the linearized mass matrices and 
gravity force vectors are written as follows: 

Mil ~ m \pl\ + m 2p^2 T ^2ph{h + 2/ 2 C 2 ) 

M P 2 = rn 2p ll 4 - miplihCi 
M p 3 = miplitj ) lie 4- TTl2p(ll + hCiWlle 
M P 4 = m lp /i^i2e 4* TKlpih + hC 2 )^\2e 
Mil ~ m 2p(h 4- h(-'2)‘ t p21e 
M P g = Tn2p(l 2 4" l\C 2 )' t P22e 

M P 2 = m2 P / 2 

M 23 = m 2 p/ 2 V’ii< ! C' 2 
M£, = m 2 p/ 2 0i 2e C 2 
M*. = m 2 p/ 2 0 2 i e 

M 26 = m 2p / 2 t/) 22e (2.46) 

= (m lp 4 - m 2 p )V>?i e 
M£, =: ( m lp + m 2 p)^ lleV’ne 

M 35 = m 2p ^ne^21eC2 

Mlz = rn2p'Pnt'l>22eC 2 
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Ml = (m lp + m 2v )ib\ 2e 

Ml 5 = m 2 pl/’l2eV’21eC' 2 
Xl P 6 = m 2 pt/’l2 e ^22eC' 2 

Mis = m 2p^L 


M 5 P 6 = m 2 pTp2U^22c 


Ml = m 2 pV' 2 2e 


G p = JTiip/iCi + na 2p (/jCi + I2C12) 

G\ = m 2p l 2 Ci2 (2.47) 

Comparing (2.46 - 47) with (2.32 - 33), six inertia parameters of mass matrices and 
gravity vectors have an analogy as shown in Table 2.1. 

Therefore, when payloads are added, six inertia parameters are changed as fol- 


lows: 

m, — ► m, 4- m, p (2.48) 

TT1 X /,£ 772|/| C 4 TTlipli (2*49) 

Di - A + m,p/, 2 (2.50) 

- LMij + m.pV’oe (2.51) 

.4M,j -► AM.j + m^Wne (2.52) 


(2.53) 
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Table 2.1: Analogy of six inertia parameters 


without payload 

with payload 

m, 

m xp 

m, :l ic 

m.p/, 

D, 

m, p lf 

LMij 

m.pi’,]' 

AMn 


NMij 



2.4.3 Centrifugal and Coriolis Force 

The velocity coupling matrix can be derived from the mass matrix using the 
Christoffel symbol (2.4). 

c _ 1 / dMi > + dM,k _ 2Mi!L\ (2.54) 

c '‘ (,) -2 { ^r + a ?J a qi 1 11 

Cjjt(i) characterizes the effects on link i which are caused by the coupled velocities 
of link j and k. The diagonal elements for j = k are the coefficients of the centrifugal 
force. The off diagonal elements for j ^ k are the coefficients of the Coriolis force. 
In equation (2.1), the states can be partitioned into the rigid body state 6 and the 
flexible body state £. 

£ Ajj + £ B,,(, + ££ PAW A + £ £ QAWA 

j= 3 j=l k = 1 J=1 *=3 


i=l 

6 6 

+ ^ ^ Rjk(i)€j£k + Gi = T t - (i — 1, 2) 

;=3 Jt=3 

£ sj, + £ D.,i +££ paw A +££ qawa 

,= 1 j=3 j=lfc=l j=lfe=3 


(2.55) 
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+ EE^w«> + t^ = ° < i = 3 ' 6 > (X56) 

j=3fc=3 J =1 


Therefore, the partitioned velocity coupling matrices can be written as follows: 

r> / -\ 1 r dAij dAik dAjk-, 

p * )m i <*, + a,, ~ a ,, 1 

( 2-57) 

n r , l.dAij dB* dB jk } 

<5 ' l(,)= 2 ( l^ + dq, dq, ’ 

(2.58) 

„ ,. x I'dBij , dB lk dD jk 
+ eq, - dqi 1 

(2.59) 

A I'dBij, dB lk dA jk 

^‘*> = 2 { aq> + aq, - 9,. ] 

(2.60) 

. 1 ,dB„ dD, k dB, k 

Q]k{ ) “ 2 dq k dq, dqi } 

(2-61) 

A 1 ,dDij , dD xk dDj k , 

aq, - a qi > 

(2.62) 


Because mass submatrix D h are not the function of the flexible body state, the 
terms related to in <?;*(*) and R jk (i ) are eliminated. The number of independent 
elements of velocity coupling matrices also can be reduced using the symmetry, the 
non-interacting, and the reflective coupling properties [<1,72]. 

C ]k (i) = C ki (i) ( 2 - 63 ) 

Cji{i) = 0 for j < i (2-64) 

Cjk(i) — —Cji(k) for j < i, k (2.65) 

However, the reflective coupling property that Tourassis and Neuman find is not 
always valid in the flexible case. Therefore, even though the symbolic manipulation 
program can be used as the computational tool, the simplication procedure must 
be accomplished under the supervision of the analyst. The symbolic programs are 
described in Appendix D. 
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Using those properties, the following independent elements of the velocity coupling 
matrix are derived. 


d\ — f {(ui e C 2 — hSi ) x 2 — (vieS? + /iC 2 )u 2 p 2 .A 2 <ix 2 } 

Jo 

d 2 i= J' ^hqnPxMdxx + 0ii e [J o (5 2 x 2 + C 2 u 2 + u le )p 2 A 2 dx 2 ] 
d 22 = J 0? 2 9«Pi-4idxi 4- 0i 2e [(^ (•S’ 2 x 2 4- C 2 u 2 4- uj e )p 2 A 2 dx 2 ] 


dz\ = 1p\u / (^ 2^2 + C2U2)p2Mdx 2 
Jo 

fi 2 

d$2 = 0l2e / (*^2^2 + C , 2U 2 P2^2^2) 
Vo 


<ill = [ 02l(02l921 +^22922)^2-^2^2 
Vo 

<f 42 = / 022(021521 + 0 22 5 22 )p 2 -4. 2 dx 2 

Jo 


2 

d 51 =2 / 021 {(“021921 + 022922) + (UleC?2 — /l 52) }p2-^2^ x 2 

Vo 

^52 ^ [ 022 {(021 921 + 022922) + (ui e C 2 - ^1 52)}p2^2^2 
Jo 

rh 

d e 1 = 0iic5 2 / ^2\p2Mdx 2 

Jo 

rh 

d$ 2 = 011e52 / tp22p2Mdx 2 
Jo 


dn = 012e52 / 021^2-^2^2 
Vo 

d? 2 = 012c52 / 022^2-^2^2 
Jo 


(2.66) 
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After evaluating the integral, 

di = (i i\ e Ci — fi5 2 )m 2 / 2c - {u\ e Si + hC2)(LM2ifa\ + LMnivi) 
d 2 \ = 0ne[ m 2^2c-?2 + (iA/ 21^21 + LM22&2)C2 + Ule m 2] + NMn^n 

d 2 2 = 0i 2 e [m 2 / 2 c 5 2 + (LM 2 161 + i -^ 2262)^2 + Ui e m 2 ] + NMuZt 2 

dz\ — '<Pue[ m 2hcS2 + (£M 2 i£ 2 1 + LA/ 22 if 22 )C 2 ] 
dj2 = V’l2e[ m 2^2c‘S'2 + (£M 2 if 2 l + -LM 22 £ 22 )C2] 

<£41 = N M2 1^21 + ( u leC2 — hSl)LM 21 
d^2 — NM22&2 "H (Ui e C 2 /lS 2 )£A / 22 

dn = NM 21&1 
d$2 = NM 22&2 

del = V>n e S 2 LM 2 i 
d 62 = 1 P 1 US 2 LM 22 


( 2 . 67 ) 


cf 71 =r \l>\ 2 eSiLM 2 l 
d>72 = ^\2eS2^^22 
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Using these elements, the velocity coupling matrices can be simplified as follows: 


C(l) = 


0 d\ ^21 ^22 ^41 ^42 

di 0 0 d 4 1 C?42 

0 0 0 0 

0 0 0 

0 0 
0 


( 2 . 68 ) 


C{ 2)« 


C(3) = 


- d j 


—d 2i 


0 ^31 

d 32 

d si 

d 52 

0 0 

0 

d 51 

d 52 

0 

0 

-<W2 

-<W2 


0 

-dn/2 

— / 2 



0 

0 




0 

—dzi 

0 0 

—d 6i 

—dg 2 

—d 3i 

0 0 

-dei/2 

-<W2 


0 0 

0 

0 


0 

0 

0 


0 0 
0 


(2.69) 


(2.70) 
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C( 4) = 


d-22 ~d-22 

0 

0 

-dr i 

-d 72 

— d 32 

0 

0 

— dri/2 

-d 72 /2 


0 

0 

0 

0 



0 

0 

0 




0 

0 


0 


C(5) = 


-d 51 

dei d 7 1 

0 

0 

^51 

dei / 2 dn/2 

0 

0 


0 0 

0 

0 


0 

0 

0 



0 

0 


0 


C(6) = 


— ds 2 

de2 

d 72 0 

0 

—dei 

<W2 

dra/2 0 

0 


0 

0 0 

0 



0 0 

0 



0 

0 


0 


2.5 Connecting and Actuating Link 


(2.71) 


(2.72) 


(2.73) 


2.5.1 Mass Matrices and Gravity Force Vectors 

Derivation of equations of motion of the connecting and actuating link is similar 
to those of the lower and upper link. Position vectors are shown in Fig. 2.4.a and 


2.4.b. 




37 


f-j = 4* x 3 cos(0i + $ 3 )]i 4* [/i r sin$i 4* Xjsin(di 4" ^ 3 ) ]j (--*4) 

r 4 = [ l lr cosd ! 4- hrCos{e x 4- 0 3 ) + xaCos{6 x 4- 0 3 4- 0* ) - 4r 9 3 + 0«)]i 

+ [i lr «^, + hrsinid, 4- 0 3 ) 4- x 4 sin(6 x + 0 3 4- 0 4 ) 4- u 4 cos(9 x 4- 9 3 4- ^)]j (2.75) 


where 

u 4 (x 4 ,f) = 041 (^ 4 )^ 41(0 4" 042 (^ 4 )^ 42(0 
Suppose generalized coordinates are defined as 

234 = {9l,93,9 4 , ^ 41 , ^42} T 

Velocity vectors are expressed by the Jacobian matrices: 

r*3 = J 3 CI 34 


(2.76) 


(2.77) 


where 


J-, = 


r*4 = J 4234 

— llrSi — X 3 S 13 — X 3 .S 13 0 0 0 
l\ r Ci 4* X 3 C 13 X 3 C 13 0 0 0 


(2.78) 

(2.79) 

(2.80) 


J 4 = 


~hrS\ — h r S\3 — U4C134 — X 4 iSx 34 ~hrS\3 ~ U4C134 X 4 Si34 

l\rC\ 4* hrSl3 ~ U4‘S’l34 4" X 4 Ci34 l3rC\3 ~ U 4 S\M 4" X 4 Ci34 


(2.81) 


— U 4 C 134 — X 4 Sj34 — ^4l5l34 — ^42 *?i34 

— U 4 5l34 4- X 4 C 134 $41^134 f42Cl34 

Mass matrices and gravity force vectors axe derived from the Jacobian matrices. 


M? = [ J 3 J3P3-A.3di3 4 ~ f J 4 J 4 p 4 Aa< 1 x 4 

1 Jo Jo 

Gi = f 3 ^ 3 ( 2 , \\p 3 Mdx 3 4- f ^ 4 ( 2 , l]p 4 A,<fx 4 

Jo Jo 

Gf = i/j[2, 2]p 3 /i 3 dx 3 4- J q * J 4 [2,2)p 4 .44dx 4 


(2.82) 


(2.83) 
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G? = [ U J 4 [2,3}p<A 4 dx4 
Jo 

Elements of mass matrices and gravity force vectors are: 

Mm = / V 3 + ^lr + 2l XT x 2 Cz)pzMdxz 

Jo 

+ [ l \ll + llr + “l + *l 

Jo 

+ ^hrihrCz + X4C34 — U4S34) 

+ 2 / 3 r (X 4 C4 ““ Ti 4 *S' 4 ))/> 4 ^4^X4 

rl 3 

MS = / (*3 + hrX Z Cz)pzA Z dx Z 

Jo 

+ / (*3r + U 4 + X 4 

Jo 

+ ^lr(^3rC r 3 + X 4 C 34 — U 4 S 34 ) 

+ 2/3r(X4 C4 — U 4 S4))p4^A^4 


<3 


= f (x\ + U 4 + hr(l4CM — U 4 S 34 ) + hriXiC* — 

Jo 

M} 4 = f 041 (X 4 + + ^3r £ 4 ) 04 -4-4^24 

= / 042( x 4 + ^ 11-^34 + 

Jo 


m 3 = jf <3 X3P3A3&3 + jTOs, + u l + x l + 2 M *« C 4 - hrU^S 0 )p^A 4 dx. 

- f‘\ X l + + U 4 + MX 4 C 4 “ U 45 , 4 ))p 4 A l dX 4 
Jo 

rl 4 

VlZ = / <M *4 + hrC A )p*A 4 dx* 

Jo 




(2.84) 
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M. 


= f ^ 42(^4 ^ 3 r ^ 4 )p 4 ^- 4^*1 

Jo 


Mg = [ U {x] + ul)p 4 A 4 dx A 
Jo 

rl\ 

M® = y xb A \x A p A A A dx A 

Mf 5 = f ip A 2X A p A A A dx A 
Jo 


M AA = f tp] x p A A A dx A 

Jo 

fl 4 

A/® = j[ IpAlVuPAAidXj 


rl\ 

Mg = J t4* 4 2P4^4^^4 


Gf = J 3 {hrC 1 + X 3 C 13) p3A 3 dx 3 

+ / (/irCi + /3rCi3 — «4‘S'l34 + I 4Cl34)/ , 4-44^ I 4 

JO 

Gf = J h X3Ci3PoA 3 dx3 + J o \hrC 13 - u 4 5i34 + x A Ci 3A )p<A A dx A (2.85) 

Gf = J (—u A S\34 A x A C\3 A )p A A A dx A 

After the integral are defined using (2.23 ■ 28), the mass matrices and gravity vectors 
of linearized equations axe : 

f® = m 2 /i r + £3 + 2 m 3 /i r / 3 cC 3 + I >4 + n^lr + ^ 3 r ) 


Af 


+2m 4 (/i r /3rC'3 + l A chrC A + hrUc^34) 

M?n = D 3 + rri 3 li r l 3 cC 3 + m A ll r + D A 
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+ m 4 (^lr^3rC , 3 + 2/ 4c /3 r C 4 + UJlrCy^) 

Mf 3 = D4 + m 4 / 4 c( ^3r ^4 + ^lrC3 4 ) 

Mj = AA /41 + £-V/ 4 i(/3rC 4 + ^lrC*3 4 ) 

Mi S = AA/ 42 + iA/ 4 2(73rC' 4 + JlrCw) 

A^22 = ^3 4" ^4^3r 4* ^4 4" 277l 4 / 4c /3 r C* 4 

jV /23 4 4“ ^n .4 ^4c^3r ^*4 

Af^ = AM41 + LM^il^rC 4 

Ml = AM 42 + LM42I3T C 4 

M32 = D4 

= AM 41 
Af^ = AAf 42 

M* = JVM 41 

rl\ 

MR = / ^ 41 042^4 A|<ix 4 

Jo 

Aff 5 = 'VM 42 

Gf = m 3 (/i r Cl + /3cCi 3 ) + rn 4 (lirC\ + /3rCl 3 + U c C\m) 

= mjlzcCw + m4(/3rCj3 + ^4cCl34) 

Gf = 


(2.86) 


( 2 . 87 ) 
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The stiffness matrices axe: 

K„ - El, fe] 1 *:, (i = 4, j - 1,2) (2.S8) 

Jo OX x 

2.5.2 Centrifugal and Coriolis force 

The velocity coupling terms are obtained from the mass matrix using the Chnstof- 
fel symbol. These terms are also simplified by several structural properties. 

t\ = /ir[ J X^SzP^A^dx^ + J {hrS3 + X 4 S 34 + U 4 C 34 ) pA^dx 4 ] 

t 2 = / {hr( x 4$34 + U4C34) + l3r{ x 4S4 + U4C4) }/> 4 ^ 4^ x 4 

Jo 

£3 = ^ 3 r[ J {X4S4 + U±C4)p4A4dX4\ 


e 4 1 = [ 041 { — ( 041&1 + 042 ^ 42 ) + (flrS 34 “ 

Jo 

e 4 2 = / ^ 42 {— (“041^41 + 042^42) + (flrS34 “ ^3r54)}p4 ^ 4 ^X 4 (2.89) 

Jo 

651 = / 041 {—(041^41 + 042^42 ) + }/>4 ^ 4 ^X 4 

Jo 

e 52 = [ 042 { (041 £41 + 042^42) + kri>4 }/?4 A 4 dx 4 

Jo 


rl\ 

C61 = / 04l(04l£ll + ^42^42 )/>4A<<ix 4 

JO 

r U 

^62 = / 042(041 £ll + 042^42 )PA^dx^ 

Jo 


After evaluating the integral, 

d = /ir[(m 3 /3 c + m 4 /3r)53 + m 4 / 4c S 34 + (LM 4 l£u + 
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e 2 = / lr [m 4 /4 C 534 + (L.V/ 4 l£u + 

+ / 3r [m4/4c54 + (LMi 1^41 + ■^■A / /4 2 ^4 2 )C4] 

e 3 = / 3r [m 4 /4c54 + (.LA^l^ll 4- iA/42^42)C4] 

c 41 = -NMnUi + (/irS 3 4 + l*rS*)LM A1 m A (2.90) 

e 42 = — jVA/ 42 ^42 + (JlrSz* + /3r-S’4)ZiiV/4 2 JTl4 

651 — — iVA/41^41 “H 

652 = — iVjW42^42 ^3r *5*4 LAI 42 ^4 

e61 = — iVjV/41^41 
e62 — “ JV’il/42^42 

Using these coefficients, the velocity coupling matrix can be expressed as follows: 

0 — e 3 — e 2 — e 4 x — e 42 

— e 2 —t\ — e 4 i —^42 

C(1) B = —Cl -641 — e 42 

0 0 
0 
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C{ 2) b = 


C(3) b = 


0 - 

-e 3 

— &51 

—652 

0 - 

-e 3 

— €51 

®52 

- 

-e 3 

-e 5 i 

€52 



0 

0 




0 

e 3 

0 


— ^62 

C 3 

0 

^61 

— ^62 


0 

^61 

— «62 


0 

0 


C(4) s = 


C( 5) b = 


eu esi e^i 0 0 

esi e6i 0 0 

e6i 0 0 

0 0 
0 

e42 e 52 e 62 0 0 

C52 e 62 0 0 

e62 0 0 

0 0 
0 


(2.92) 


(2.93) 


(2.94) 


(2.95) 


2.6 Conclusion 

Mass matrices and gravity vectors axe directly derived from the Jacobian matri- 
ces which axe easily calculated from position vectors by SMP. Because the deriving 
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procedure is simple, it reduces the possibility of producing incorrect equations. Fur- 
thermore, this form can easily be used to expand a series of mode shape functions 
describing elastic deformation. Six inertia parameters are defined and analogous 
terms exist for the cases with and without payload. The velocity coupling matri- 
ces, which are the coefficients of centrifugal and Coriolis force terms, are derived 
from the mass matrices using the Christoffel symbol and are simplified using sev- 
eral structural properties. The resulting velocity coupling matrices have a structure 
which is useful to reduce the number of terms calculated, to check correctness, or 
to extend the model to higher order. 
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CHAPTER III 

Constrained Dynamic Analysis 

3.1 Introduction 

In this chapter, a numerical method is presented for dynamic analysis of closed 
kinematic chain systems such as parallel link mechanisms. Dynamic solution of a 
closed kinematic chain system requires solution of a mixed set of differential equa- 
tions of motion and algebraic constraint equations. Singular Value Decomposition 
(SVD) is used for reducing the equations by eliminating the unknown constraint 
force. The nonlinear constraint equations of RALF are derived. Natural frequencies 
and eigenvectors of a closed kinematic chain system are derived and verified by a 
simple example. 

3.2 Singular Value Decomposition 

As mentioned in chapter II, the equations of motion of a flexible arm with 
a parallel link mechanism are expressed by the mixed set of differential equations 
(2.6) and nonlinear algebraic equations (2.7). The Lagrange multiplier A is included 
in differential equations to describe the unknown constraint force by the constraint 


Jacobian matrix. 


Mq + Kq + F + $*\ = Q 


(3.1) 
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where F is the nonlinear force vector which includes Coriolis and centrifugal forces 
and the gravity force. Nonlinear algebraic constraint equations describe the rela- 
tionship among the angles within the closed kinematic chain. 

$(9)= 0. (3.2) 

Differentiating the constraint equation (3.2) with respect to time yields the velocity 
form of constraint equations. 

$ q q = 0 (3.3) 

For solving this mixed set of differential and algebraic equations, the unknown 
constraint force vector, A, has to be eliminated from the differential equations. The 
constraint Jacobian matrix with rank m, can be decomposed into the following 
form using Singular Value Decomposition (SVD). 

= UZV T (3.4) 

With proper partitioning [38], it can be expressed as 

V T 

<t> q = [U, V 2 ][S m 0] 1 (3.5) 

V? 

where U : and Vj are orthonormal bases for four fundamental subspaces as shown in 
Fig 3.1. The columns of £7, are the orthonormal eigenvectors of the matrix 
The columns of V are the orthonormal eigenvectors of the matrix $^$ 7 . is equal 
to diagCcTi,^,...,^) ordered <r a > <r 2 > ... > 0. The are called singular values 
which are the nonnegative square roots of the corresponding eigenvalues. Notice 
that V 2 is the null space of $ 7 which satisfies the following relationship. 

$,V 2 = 0 


(3.6) 
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Also, called the pseudo inverse of <£,, is defined as follows: 

T" 1 1 [ U? 

= [Vi V 2 ] (3-7) 

° J [ U? 

Premultiplying (3.1) by V 2 and using the orthogonality of (3.6) gives 

V?Mq + V?Kq + V 2 t F = V?Q (3.8) 

The unknown constraint forces axe thus eliminated from the equations of motion. 
However, because V 2 M is the (n - m) x n rectangular matrix, additional equations 
are needed to get a solution. Let us define a new variable z which is an independent 
coordinate with dimension (n - m) x 1. The homogeneous solution to (3.3) is V 2 z. 

q = V 2 £ (3.9) 

Geometrically, it is the projection of the velocity vector q onto the tangent hyper- 
plane of the constraint surface. Moreover, the time derivative of (3.3) gives 

$ q q = 

= -i r (* t ),ff (3-io) 


Then, q is represented as 

q = -$qq T ($q)qQ + V 2 Z (3.11) 

The first and second term in the right hand side of (3.11) are the particular solution 
and the homogeneous solution to equation (3.10) respectively. Physically, they are 
the normal accelerations and the tangential accelerations of the constraint surface, 
respectively. Finally, by integrating (3.9), position vectors are expressed as 


q — V 2 : + C 


(3.12) 
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where C, a constant, is here chosen to be zero to satisfy initial conditions. 

Substituting (3.9), (3.11), and (3.12) into (3.8), the equations of motion become: 

VfMVtz + V 2 t KV 2 z = V 2 t Q + V 2 t F (3.13) 

where 

F = M$+q T {$ q ) q q ~ F ( 3 - 14 ^> 

Equation (3.13) is a set of n - m equations in terms of the independent generalized 
coordinate z. As a result, the n equations (3.1) and m constraint equations (3.3) 
are reduced to n - m equations in (3.13) by the coordinate transformation matrix Vj 
[67]. Because the independent generalized coordinate lies on the tangential plane 
of the constraint surface instantaneously, the changes in generalized coordinates 
due to the integration of i during a small time interval do not result in significant 
constraint violation [37]. Therefore, the reduced equations of motion (3.13) are free 
from constraints and stable for numerical integration error. 

3.3 Constraint Equations 

To apply the SVD method to RALF, the constraint equations must be described 
first. Suppose the elastic deformation is small compared to the length of the link, 
the deformed parallel link mechanism is depicted as in Fig 3.3. For the virtually 
cut joint C", the constraint equations axe expressed by two vectors. 

AB + BB 1 + B~C' = AD + DC' (3.15) 


or 


( l\jC\ — Ui e Si + ^rCn)* + + U\ e C\ + 

— (^ 3^13 "b 14^134)1 + (I3S13 "h ^4^134^ 


(3.16) 
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& = 


= 0 


(3.17) 


When this equation is expressed in the form of equations (3.2) and (3.3), the position 
form of constraint equations are 

h/C\ — U\ e S\ + hrC \2 — hCu — /4C134 
h/Sl + U\ e C\ + UrS \2 — US \3 — /45134 
or, in the velocity form of constraint equations, the elements of the constraint 
Jacobian matrix are: 

c 

UjS\ + u \eC\ + hrSll — I3S13 ~ US \34 

hjCi — uj e Si + UrC 12 ~ hC \3 — UC134 

hrS \2 —hSl 3 ~ USl 34 ~ US \34 

UrC 12 —UC 13 — U Ci 34 —UC134 


«{ = 


(3.18) 


i>UeSi V’HeSl 0 0 0 0 

V’lleC'l 4 >UeC\ 0 0 0 0 

In the rigid parallel link mechanism as shown in Fig 3.2, the elastic deformation 
is not included in the constraint equations. The coordinates transformation matrix 
V 2 , which is derived from by SVD, is 

UjS\ + UtS\2 ~ US\3 ~ US\3\ 

h/Ci + UrC 12 — I3C13 — I4C134 




UrS \2 ~USl 3 ~ ^4‘S'l34 —US 134 

UrC 12 — /3C13 — I4 C134 —I4C134 


(3.19) 


V 2 = 


1 

0 


1 

0 

0 

l/v/3 


0 

0.577 

0 

l/y/3 


0 

0.577 

0 

-l/y/3 


0 

—0.577 


(3.20) 
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where 

5r = {#i, 02, 03,8 4 } t (3-21) 

V 2 is independent of angles and the ratio of the link length. The independent coor- 
dinates of the rigid case are 6 X and ^(0 2 + 0 3 - 6> 4 ) for all configurations. However, 
if am independent generalized coordinate is selected by Gaussian elimination, it may 
change frequently [37]. On the other hand, in the flexible parallel link mechanism as 
shown in Fig 3.3, V 2 depends on the elastic deflection and its configuration, unlike 
the rigid parallel link mechanism, as shown in equations (3.22 - 3.24). 

For u Xe /l\f = 0, 0 2 = 135 deg. 

_1 0 0 00000 

0 0.583 0.169 -0.010 0 0 0 0 

0 0.583 0.169 -0.010 0 0 0 0 

0 -0.566 0.388 -0.023 0 0 0 0 

0 0.026 0.890 0.006 0 0 0 0 

^ 2 = 

0 0.002 0.006 0.999 0 0 0 0 

0 0 0 00000 

0 0 0 00000 

o 0 0 00000 

o 0 0 00000 
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For ui e /l\f = 0.01, 


V 2 = 


For uu/lif = 0.0, 


V 2 = 


where 


d 2 = 135 deg. 

-1 0 0 

0 0.5S0 0.184 

0 0.580 0.173 

0 -0.572 0.381 

0 0.012 0.890 

0 0.001 0.006 

0 0 0 

0 0 0 

0 0 0 

0 0 0 

9 2 = 60 deg. 

-1 0 0 

0 0.584 0.163 

0 0.584 0.163 

0 -0.572 0.394 

0 0.035 0.890 

0 0.002 0.006 

0 0 0 

0 0 0 

0 0 0 

0 0 0 


0 0 0 0 0 

-0.011 0 0 0 0 

- 0.011 0 0 0 0 

- 0.022 0 0 0 0 

0.006 0 0 0 0 

0.999 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

- 0.010 0 0 0 0 

- 0.010 0 0 0 0 

-0.023 0 0 0 0 

0.006 0 0 0 0 

0.999 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 


1 , = {«., 02, 03, 04>£n.£l3,f31> £22, £■»!,£«/ 


(3.23) 


(3.24) 


( 3 . 25 ) 
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The angle 0 4 is mainly changed by the deformation of the lower link. Because the 
links are vibrating during motion. V 2 is also changing. Therefore, the constraint 
Jacobian matrix must be decomposed at every step to preserve the direction conti- 
nuity of the basis of V 2 . However, performing S\'D at each integration step would 
result in a tremendous increase in computing time. Furthermore, performing rede- 
composition and recomputing the initial conditions on z at each time step requires 
numerical integration algorithms which use information only from the current time 
step even though this eliminates the need for null space continuity. Therefore. Mam 
introduced the velocity norm as the criterion for redecomposition [46]. V 2 is held 
constant until the specified criteria Eire violated. In order to circumvent the rede- 
composition, updating algorithm using QR decomposition was developed [37]. For 
more accurate cind stable solution, several algorithms have been developed recently 
[39,32,57]. But, those methods are not applied in this thesis. 


3.4 Computational Algorithm 


An algorithm for solving the equation (3.13) is summarized as follows: 

1) The constraint Jacobian matrix $<, is decomposed by the singular value de- 
composition subroutine LSVDF of IMSL [33]. Then, the initial condition for inde- 
pendent coordinate 2 , at time step i is defined by the following transformation. 



(3.26) 

(3.27) 


Initicd conditions for first-order differential equations are 

Qo 

-0 
z o 


Y(t 0 ) = 


(3.28) 
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2) Integrate (g,, i,-, i\) of equations (3.9) and (3.13) simultaneously from f,- to 
t, +1 to get (< 7 i + i, r, +1 , ~ 1+1 ) using the subroutine DGEAR of IMSL which is based 
on Adams predictor - corrector method. 

3) Solve the original velocity vector <],+\ from i;+i using the transformation of 
equation (3.9). 

4) Check the velocity norm. If || i || 2 is less than a predetermined fraction of 
|| q || 2 , the constraint Jacobian matrix needs to be redecomposed [37,46] and 
step 1) is repeated. Otherwise, step 2) is executed next. 

5) repeat the above steps until the final time is reached. 

Application of this algorithm to RALF will be explained in Chapter V. 

3.5 Natural Frequencies and Mode Shapes 

Natural frequencies and mode shapes can be derived from the linearized equa- 
tions of motion. Eigenvalues and eigenvectors of the constrained equations can be 
obtained from the reduced equations (3.13) linearized about zero velocity. 

V 2 t MV 2 z + V? KV 2 z = 0 (3.29) 

Eigenvalues of the constrained equations, (3.1) and (3.3), axe the same as those of 
the reduced equations (3.13). Eigenvectors of the constrained equations axe derived 
by transforming those of the reduced equations as follows. 

q = V 2 z (3.30) 

The validity of the above theory is demonstrated by the following simple example. 

Now, the flexible parallel link mechanism is somewhat analogous to two masses 
connected to one another by springs as shown in Fig 3.4. This system can be fully 
described by the two coordinates xi and x 2 . Equations of motion of this system 
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Figure 3.4: Constrained system 
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Figure 3.5: Unconstrained system 
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are: 


Mi ■+• M] M 2 
M 2 Mi 

If the values of mass and spring axe arbitrarily assumed as 
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the equations of motion are 

3 2 

9 9 

L “ “ 

and 
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A\ = l 
K 2 = 2 , 


1 0 
0 2 


= 0 


M; x K a = 


1 -2 
-1 3 


The characteristic equation is 


A 2 - 4A + 1 = 0 


r x = 


(3.31) 


(3.32) 


(3.33) 


(3.34) 


(3.35) 


(3.36) 


and the natural modes are 

r 1 -0.7321 

0.3660 1 

The system can be also expressed by three coordinates Xi, x 2 , x 3 and a constraint 
equation as shown in Fig 3.5. The redundant coordinate x 3 describes the relative 
motion of the two masses. Equations of motion in the new coordinates system, q = 
{xi, x 2 , x 3 }, are 


+ $JA = 0 
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(3.37) 
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and the constraint equation is 

x 3 = x x + d 

The velocity form of constraint equation is 

i 3 = i\ 

and it can be expressed as follows. 


/ 

i 1 


= 


1 0 -1 


x 2 

*3 


= 0 


The constraint Jacobian matrix, is decomposed by SVD. 

l/y/2 -l/y/2 0 

0 0 1 

-l/y/2 -l/y/2 0 


-l/y/2 0 
V 2 = 0 1 

—l/y/2 0 

The orthogonality of and V 2 is checked as follows : 

— l/y/2 0 



y/2 0 0 



0 0 0 



0 0 0 


Here, the null space, V 2 , is 

r 


* q V 2 = 


1 0 -1 


= [0 0 ] 


(3.38) 


(3.39) 


(3.40) 


(3.41) 


(3.42) 


(3.43) 


0 1 

-l/y/2 0 

Then, the reduced mass and stiffness matrix are used to find the eigenvalues and 
the eigenvectors. 


{V? MV 2 )~\Vj KV 2 ) = 


1 2v/2 

l/y/2 3 


(3.44) 
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The characteristic equation is 


A 2 - 4A + 1 = 0 


(3.45) 


and the eigenvectors are 


I\ = 


1 1.0353 

-0.2588 1 


The eigenvectors of the constrained system can be derived using V 2 . 


r, = v 2 r z = 


-1/^2 

-l/v/2 


1 

-0.7321 

-0.2588 

0.9659 

= 

0.3660 

1 

— 1/ n/2 

-1/V2 ' 


1 

-0.7321 


The first and second row are the eigenvector of ii and x 2 


r,= 


1 -0.7321 

0.3660 1 


(3.46) 


(3.47) 


(3.48) 


Therefore, the characteristic equation (3.35) and eigenvectors (3.36) of the second 
order system are the same as those equations (3.45) and (3.48) of the reduced system 
respectively. 


3.6 Conclusion 

Parallel link mechanisms with rigid links have a simple relationship between the 
angles of the closed kinematic chain. However, the parallel link mechanism with 
flexible links requires nonlinear constraint equations to describe the relationship 
between angles because the link deflection gives a perturbation to the relationship. 
Therefore, a mixed set of differential and algebraic equations must be solved si- 
multaneously. The transformation matrix V 2 , which is derived from the constraint 
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Jacobian matrix by SVD, is used to obtain the reduced set of equations of motion 
which have no additional constraint equations. The reduced equations have been 
solved without any significant constraint violation. Therefore, SVD is a stable and 
efficient numerical method for the closed kinematic chain system. The simple ex- 
ample shows that natural frequencies and eigenvectors of the reduced equations are 
the same as those of the constrained equations. 
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CHAPTER IV 


Mode Shape Functions of RALF 


4.1 Introduction 

Mode shape functions are employed to describe the elastic deformation of the 
flexible manipulator. In the assumed mode method, mode shape functions need 
only to be admissible functions which satisfy the geometric boundary conditions 
and form a basis set. However, a large number of modes are required to obtain ac- 
curate frequencies. The number of modes to be included can be reduced by choosing 
appropriate functions which satisfy static equilibrium at the interface between links. 
These functions can be derived systematically by using the component mode synthe- 
sis. Various component mode synthesis approaches have been developed depending 
on the assumed boundary conditions. Different boundary condition assumptions 
require different coordinate systems to describe the elastic deformation. Therefore, 
the method which fits with the current coordinate system is chosen. In order to 
explain component mode synthesis, an L - shaped beam is analyzed first. Secondly, 
mode shape functions of RALF are derived. Natural frequencies and mode shape 
functions of the analytical model of RALF are compared to the results of a finite 
element method, component mode synthesis and experiment results for validation 
of the proposed mode shape functions. 
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4.2 Component Mode Synthesis 

A complete structure can be regarded as an assembly of component structures. 
The dynamic behavior of each component can be represented by a reduced number 
of component modes. Modes of the complete structure are formed from the compo- 
nent modes by enforcing equilibrium and compatibility along component interfaces. 
Generally, the constraint modes are introduced to provide shear and moment com- 
patibility when the structure is assembled. Therefore, component mode synthesis is 
useful for predicting accurately the lower modes and frequencies of a structure us- 
ing a small number of component modes. Numerous methods for component mode 
synthesis have been presented. Fixed, free, or loaded boundary conditions are used 
to determine the component modes. An L - shaped beam as shown in Fig. 4.2.a is 
used to explain the difference between these methods. 

In a fixed interface method [19], modes of each component are expressed by a 
constraint mode and a normal mode. The constraint mode is defined as a static 
mode of internal nodes resulting from an imposed unit displacement at the interface 
node. The normal mode is defined as a dynamic mode of internal nodes when the 
boundary condition is in effect on the external nodes. In Fig. 4.1. a, link deflections 
of the lower beam can be expressed by rigid body rotations and clamped - clamped 
modes. Here the rigid body rotation is the constraint mode and the clamped - 
clamped modes axe the normal modes. The boundary conditions of the lower beam 
and the upper beam are clamped - clamped and clamped - free, respectively. 

In a free interface method [27], modes of each component are expressed by 
normal modes only. Modes of the upper beam are defined from the mass center. In 
Fig.4.1.b, the boundary conditions of each beam are clamped - free and free - free. 

In a loaded interface method [6], modes of the lower beam sure modified by 



63 




Figure 4.1. a: Coordinate system of a fixed interface method 



Figure 4.1.b: Coordinate system of a free interface method 



Figure 4.1.c: Coordinate system of a loaded interface method 
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added inertia and stiffness which are the effects of the upper beam. In Fig. 4.1.c, 
the boundary conditions of each beam are clamped - mass - spring and clamped - 

free. 

Even though all methods give the same result in the limiting case, the coordi- 
nates which are used to describe the elastic deformation are different. Coordinates 
of the loaded interface method match closely with the coordinate system of RALF 
as seen in Fig. 2.3. The x, axis is aligned with the tangent of the respective link at 
the origin 0 X . Therefore, the mode shape functions of each link are derived using 
the loaded interface method. 

4.3 Loaded Interface Method 

In this section, a loaded interface method [6] is summarized. An L-shaped beam 
as shown Fig. 4. 2. a is an example configuration of a two serial link manipulator and 
the typical example of component mode synthesis [31]. In this section, component 
mode synthesis is explained using the L - shaped beam. In Fig. 4.2.b, the complete 
structure can be separated into two components - the main component and the 
branch component. Each component’s modes are expressed in terms of internal (i) 
coordinates and junction (j) coordinates. 

A constraint mode is defined by statically imposing a unit displacement at the 
junction node and zero displacement at the internal node [19] - 


Ku 

Ka 

< 

4>'b 


0 


K» . 


I 

* 


R 


where R is the reaction at the junction. Thus, the constraint mode is obtained from 
the top row partition. 


(4.2) 
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A normal mode is obtained from the internal coordinates only. 


(K u -u\\l x ,w£ =0 


(4.3) 


The normal mode set is truncated to a set of normal modes <f>' B . Therefore, 
the internal node displacements of the branch component B can be expressed by 
superpositions of constraint modes <j>' B and normal modes 4>g • 


u' B = 4'b u b + &'b £b 

The total node displacements of branch component B are 


(4.4) 


UB = < 


U, 


u\ 
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> = 
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0? 4>\ 



(4.5) 


The total potential energy for uncoupled components A and B can be written 


as 


1 T 1 T 

PE = -u T A K A u A + -u b K b u b 


(4.6) 


Components A and B axe coupled together by constraining the interface coordinates. 
This constraint is expressed as 

(4-7) 


u'b 


Equation (4.7) assumes that the same reference coordinate system is used for both 
components. Different coordinate reference systems will require a rotational coor- 
dinate transformation using directional cosines. Therefore, the coordinate transfor- 
mation that couples components A and B together becomes 
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(4.8) 


Therefore, equation (4.6) can be rewritten as 

PE = \u t a K a u a ( 4 -9) 

2 

where 

K a = K a + T?K b T, ( 4 - 10 ) 

In similar way, the kinetic energy is 

KE=\uImK a (“-id 

where 

M a -Ma + Tj M bT\ ( 4 - 12 ) 

As shown in Fig. 4.2.c, the mass matrix M A and the stiffness matrix Ka of compo- 
nent A are modified by the static influence of the branch component B. 

Using these modified matrices, equations of motion of main component A are 

M a u a + K a u a = 0 ( 4 -13) 

The eigenvector of (4.13) yields 

u A = ( 4 - 44 ) 


Eigenvectors *1 * a are also truncated to a set of normal modes O a • 

u 


u. 


<4 


*'a 
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U A = \ 


(4.15) 
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Both components are coupled using compatibility of junctions. 

U B = U A = 4aZa ( 4 ‘ 16 ) 

So, 

k ] I 0 

Ug * = 4>{\ 0 

. 0 / 

Therefore, using equations (4.15), (4.5), and (4.17), node displacements of each 
component can be expressed in terms of reduced order modal coordinates. 





(4.18) 


The reduced order of the mass and stiffness matrices axe derived by the trans- 


formation matrix T 3 . 
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T. 

M b 
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T. 

K b 


(4.19) 

(4.20) 
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Natural frequencies and modes of the complete system are calculated from eigen 
values and eigenvectors of the following equations. 



The node displacements of each component are derived from eigenvectors of Eq. 
(4.21) using transformation matrix T s of Eq. (4. IS). 

4.4 Component Modes of RALF 

For the validity of modes which are derived using component mode synthe- 
sis, system natural frequencies of finite element model, component mode synthesis 

model, and the analytical model are compared. 

The finite element model of RALF is modeled as shown in Fig. 4.3. a. Note that 
for this model the lower link is assumed to be supported where the connecting link 
attaches. This is a simplification over the physical system. The support sleeves 
of the lower link and the upper link and the connecting link are assumed as rigid 
elements by setting Young’s modulus of these elements 1000 times bigger than that 
of the other elements. Because the hydraulic actuators are also assumed rigid, the 
boundary conditions of node A and C are assumed to be simply supported as well 
as node O.'The finite element model is simulated by MSC/PAL2 program [51]. 

Component mode synthesis can predict the lower modes and natural frequencies 
of RALF using a small number of component modes. The derivation procedure is 
similar to Sec. 4.3. However, a little modification is required because RALF has a 
parallel link mechanism with components connected by pins. 

RALF can be partitioned into three components as shown in Fig. 4.3.b. First, 
the boundary condition at node D' of the actuator link is assumed as simply sup- 






ported imposing a fixed boundary condition on the pin joint. The reaction force 
or torque can be expressed by adding the stiffness and the mass at node D of the 
upper link using (4.10) and (4.12). Next, the boundary condition at node B of the 
upper link is also assumed as simply supported. Furthermore, the stiffness and the 
mass are also added at the node B of the lower link. Using these boundary condi- 
tions as shown in Fig. 4.3.C, the node displacements of each lmx can be obtained b\ 
the finite element method. Node displacements of each component of the uncoupled 
system are expressed as follows. 
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The upper link has two interface nodes with the actuating link and the lower link. 

uf = ui = KU ( 4 - 23 ) 

u J c = u [ 2 = 4 >i 2N £b + 4>i c u[ x 
= 4>l 2N t h + 4>l c 4>ii« 

= + V’flso (4.24) 


Using (4.23) and (4.24), uncoupled generalized coordinates can be described in 
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terms of generalized modal coordinates. 
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(4.25) 


Relationships between node displacements and generalized modal coordinates are 
derived from (4.22) and (4.25). 
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(4.26) 


In the pin joint connection, equations (4.23) and (4.24) are modified as follows. 
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Jl = ui = = <t>ii« 


(4.27) 


u 3 c = uf = T p (0f 6 + 4>l C T p ul l ) 
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because each node is defined in three coordinates - 2 translational and 1 rotational. 
Therefore, equation (4.26) is also modified as follows. 
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Using this coordinate transformation matrix T,, the reduced order of mass and 
stiffness matrices are derived using (4.19) and (4.20). Natural frequencies and the 
node displacements of RALF can be derived from this reduced order mass matrix 
and stiffness matrix. 

In chapter 2, the following modal data of each component are required for the 
equations of motion. 

NM =s j ip 2 dm 

= ( 4 . 31 ) 

LM = J dm 
AM = J xxfrdm 

where dm is the mass of the small segment. These data can be obtained by several 
methods. 

First, if the values of the interface mass and spring are known, these data can be 
obtained from the analytical solution. However, it is not easy and it is a time 
consuming job to get an analytical solution. 

In a second method, the mode shape functions can be derived from node displace- 
ments of the finite element model using a polynomial. However, this method has 
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a problem in calculating the stiffness matrix because the second derivative of the 
estimated polynomial has some difference from the true value. That is, it is not 
practical to find a polynomial to satisfy the displacement and slope of every node 
point at the same time. 

Bv the third method used here, modal data can be expressed by summation instead 
of integration. Therefore, modal data can be obtained from node displacements of 
lumped mass model of each component. 

NM = 

KW = * t K* (4-32) 

jv 

LM = ^2 TKjVii 1 j) 

1=1 

N 

AM = y 2 x jmjibi( x j) 

1=1 

where 

M = unloaded mass matrix 

K = stiffness matrix 

<5 = mode shape matrix 

mj = mass of the j tk finite element node 

=i mode shape of the i tk mode at the j th node 
x ■ = vector locating j th node with respect to the local frame 
Mode shape matrices of the lumped mass model are less reliable than those of 
the consistent mass model. Therefore, it was observed during a simulation that 
the mass matrix is not positive definite for a model with small number of nodes. 
However, modal data can be readily and inexpensively computed from the lumped 
mass model which has a large enough number of nodes. 
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Table 4.1: Comparison of Natural frequencies 


Modes No. 

FEM 

CMS3 

CMS2 

AM3 

AM2 

AMI 

1 

S.3S 

S.39 

S.43 

8.41 

S.45 

S.46 

2 

15.40 

15.55 

15.72 

15.61 

15.98 

16.20 

3 

30.51 

30.52 

30.54 

30.53 

30.54 

30.56 

4 

92.49 

94.10 

101. SI 

9S.62 

106. 5S 

10S.61 

5 

117.60 

119.06 

123.67 

119.88 

121.S5 

122.65 

6 

120. SO 

120. S4 

121.12 

120.S9 

121.22 

121.33 


The first six natural frequencies of three models - a finite element model, a 
component mode synthesis model, and an analytical model - are compared as shown 
in Table 4.1 . In this table, component mode sets of the component mode synthesis 
(CMS) model and the analytical model (AM) are truncated at two mode (CMS2, 
AM2) and three mode (CMS3, AM3) per link. If three modes are included, all six 
natural frequencies match each other quite well. On the other hand, if two modes 
are included, only the first three natural frequencies match each other to within 
10 %. However, because we axe interested in the lower frequencies, two modes are 
enough for the reduced order analytical model. Furthermore, AMI is the model 
which uses clamped - mass, clamped -free, and pin - pin boundary conditions for 
the lower, the upper, and the actuator link. This model is more different from 
the FE model than AM2. Therefore, component mode synthesis is the systematic 
method which can derive proper mode shape function. 

4.5 Discussion of Experiment 

In fact, there are discrepancies between the analytical model and the real system 
because the real system has a complex structure which is difficult to analyze by the 
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Table 4.2: Comparison of natural frequencies of the modified model 


Modes No. 

FEM 

CMS3 

CMS2 

AM3 

AM2 

AMI 

1 

6.03 

6.04 

6. OS 

6.06 

6.09 

6.11 

2 

15.24 

15.30 

15.49 

15.37 

15.74 

15.99 

3 

30.74 

30.74 

30.75 

30.75 

30.76 

30.79 

4 

75.63 

77.21 

S4.91 

81.73 

89.70 

91.73 

5 

98.25 

99.70 

101.3 

100.5 

102.5 

103.3 

6 

120.2S 

120.32 

120.59 

120.37 

120.70 

120. SI 


analytical method. The sectional areas of beams axe not uniform. The links are 
connected with offset brackets. The connecting joint axis between the hydraulic 
actuator and the lower link is not colinear with the axis joining the connecting 
link to the lower link as shown in Fig. C.l of Appendix C. Therefore, the flexible 
part of the lower link increases. Table 4.2 shows the natural frequencies of each 
model of the modified structure. The finite element models are used to explain the 
discrepancies between the analytical model and the real system. Two types of finite 
element models are created. One is a simplified model using beam elements with 
seven different El (modulus of elasticity times axea moment of inertia) values and 
one lumped mass element at the end of the lower link. The simplified model assumes 
support sleeves and connecting link as rigid similar to those used in the analytical 
model. The other finite element model is a detailed model using beam elements in 
which thirteen different El values and three kinds of lumped mass elements were 
used. The detailed model used makes no assumptions about the rigidity of any of 
the links. In experiments, an electromechanical shaker was attached to the structure 
at the end point of the lower link (Point B in Fig. 4. 3. a). Because turnbuckles were 
used in place of the hydraulic cylinders, the boundary condition can be assumed 



fixed. The detailed procedures of the experiment are described by Huggins [28]. 

Figure 4.4 shows the mode shapes of RALF and the associated system natural 
frequencies. The first system mode is dominated by the bending of the lower link. 
The frequencies of the first mode are nearly equal each other. The second system 
mode is dominated by the bending of the upper link. However, there is a discrep- 
ancy between the frequencies of the second mode of the simplified model and the 
detailed model. The boundary condition of the upper link is more complicated in 
the real system. Therefore, a more accurate model is required to predict the true 
frequency. The third system mode is dominated by the pinned - pinned bending 
mode of the actuator link. However, there is a discrepancy between the predicted 
and the measured natural frequencies. A large amount of friction in the pin joint 
causes the joint to exhibit some characteristics of a clamped end condition. This 
frequency is measured when the structure was excited by a shaker and the joints are 
fixed by a turnbuckle with no bearing. However, if the structure is excited by the 
hydraulic cylinder, the frequency of the actuator link decreases. Therefore, the rea- 
son for the frequency difference in the third mode appears to be friction at the joint 
which is reduced by actuator motion. The higher frequencies and mode shapes have 
discrepancies between the analytical models and the finite element models because 
the analytical model included only two modes per link for describing link deflec- 
tions. The higher modes and frequencies of the experiment also have discrepancies 
because of their small signal to noise ratio, out of plane motion, and the closeness 
of the frequencies as shown in Fig. 4.5. Even in the finite element model, the higher 
modes axe sensitive to modeling. 
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Figure 4.4: Comparison of mode shapes and natural frequencies 

between analytical model, finite element model 
and experiment 
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Figure 4.5: Frequency spectrum of RALF 
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4.6 Conclusion 

In the assumed mode method, there are numerous possible choices of admissible 
functions. However, a laxge number of modes are required for the completeness of 
the set of admissible functions. In order to speed up the convergence, mode shape 
functions which satisfy the static equilibrium are used in component mode synthesis. 
Among several component mode synthesis approaches, the loaded interface method 
fits with the current coordinates system. Proper mode shape functions for a reduced 
order model are obtained using the loaded interface component mode synthesis. 
Comparison between the finite element model, the component mode synthesis model 
and the analytical model show that mode shape functions which are determined by 
component mode synthesis improve the convergence. Because the real system has 
a complex structure which does not match certain simplification necessary in the 
analytical model, there are some discrepancies between the analytic model and the 
experiment. These discrepancies are explained using a simplified and a detailed 


finite element model. 



82 


CHAPTER V 


Verification of Analytical Model 


5.1 Introduction 

The experimental verification of nonlinear dynamics is a difficult job which has 
not been extensively studied compared to linear dynamics. Possible verification 
methods include a time domain and a frequency domain methods. However, the 
time response with nonlinear dynamics depends on initial conditions and input am- 
plitudes. Therefore, an accurate dynamic model is required. As mentioned in the 
previous chapter, the real system has unmodelled dynamics and parameter uncer- 
tainties hard to describe analytically. Hence, direct comparison between the analyt- 
ical model and the real system was difficult. In this chapter, several techniques are 
tried to verify the nonlinear dynamics of RALF. First, as an alternative method for 
a time domain verification, a TREETOPS model is used for the verification of the 
analytical model. Secondly, as one of the frequency domain methods, the nonlinear 
dynamics of RALF is studied by a sinusoidal excitation. Finally, actuator dynamics 
has a significant interaction with flexible body dynamics. The effect of hydraulic 
cylinder dynamics on flexible body dynamics is discussed. 
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5.2 Analytical Model Verification using TREETOPS 

TREETOPS is a computer simulator of the dynamics of a flexible multibody 
structure with loop closures. TREETOPS is being developed by DYNACS Inc. 
under contract to NASA’s Marshall Space Flight Center (MSFC). The program 
was obtained from Dr. Henry B. Waites of MSFC. The name TREETOPS refers to 
the class of structures which may be simulated by the program, that is those having 
a tree topology. A tree topology is one in which multiple bodies are connected 
by rotational and translational joints. TREETOPS can also handle the structure 
with loop closures by constraining the position or velocity of connecting joints. The 
primary output of the program is a time history of flexible body motion in response 
to an active control system consisting of actuators, sensors, and controllers. The 
equations of motion are derived using Kane’s method. The TREETOPS model is 
described in Appendix E. 

The step responses of the analytical model and the TREETOPS model are 
compared as shown in Fig. 5.1. The time responses of the two models match each 
other fairly well even though they use different methods for deriving the equations 
of motion. The constraint error measured at the connection joint is stable as shown 
in Fig. 5.2. The analytical model used the SVD algorithm to generate the matrix V 2 
while TREETOPS model used the QR algorithm. SVD algorithm updates the V 2 
matrix based on the given criterion. On the other hand, the QR algorithm updates 
the V 2 matrix at every time step [37]. Even though the constraint error history 
using the SVD algorithm is a little different from the constraint error history using 
the QR algorithm as shown in Fig. 5.2, time responses of other variables are not 
different as shown in Fig. 5.1 . The constraint error also depends on the integration 
method. Because the flexible body dynamics produces a stiff system which has 


84 


several widely spaced frequencies in one system, Adam s method is used. The 
Runge-Kutta method allows the constraint error to increase as shown in Fig. 5.3. 
The simulation time of analytical model is shorter than that of the TREETOPS 
model. For example, the CPU time of analytical model is about 6 minutes in 
VAX/750 for 1 second time response and 1 millisecond time step. On the other hand, 
the CPU time of the TREETOPS model is about 35 minutes. The expected reason 
for this difference is that TREETOPS is a general program which can handle a 3 
dimensional motion of a multi body. Furthermore, the dynamic equations of motion 
are calculated recursively by matrix manipulation. Whereas, the analytical model 
is expressed in symbolic form explicitly. In both models, the system frequencies 
which are related to the upper link are high compared to the experiment results. 
This difference comes from the actuator dynamics effects on the system dynamics. 
The hydraulic actuator dynamics will be discussed in the section 5.4. 

5.3 Verification of Nonlinear Dynamics 

5.3.1 Analysis of Nonlinear Dynamics 

In order to verify the nonlinearity, the following questions are expected. Where 
does the nonlinearity come from ? Which excitation method is adequate to test the 
nonUnear dynamics ? Which sensor is adequate to measure the nonlinear effect ? 

First, nonlinear dynamics of a flexible manipulator comes from configuration 
change, Coriolis and centrifugal force, and the coupling between rigid body motion 
and flexible body motion as shown in the equations of motion. Secondly, a step time 
response has been the common method for model verification. However, even though 
the step input can excite many modes at one time, the time response is susceptible 
to the unknown dynamics, such as the signal noise and the friction. Furthermore, it 
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is difficult to detect the amount of nonlinearity from the step response. On the other 
hand, the frequency spectrum of the response to a sinusoidal excitation is a method 
which can detect the nonlinearity. The response of a nonlinear system is dependent 
on the amplitude of the excitation. One advantage of a sinusoidal excitation is that 
the input force can be precisely controlled. Furthermore, large amounts of energy 
can be input to the structure at each particular frequency. Therefore, it results in 
relatively high signal-to- noise ratios and allows the study of structural nonlinearities 
at any specific frequency. Third, the time response can be measured by several 
sensors. The possible sensor types axe accelerometer, encoder, tachometer, and 
strain gauge. 

Before the experimental verification, the simulation result can show the general 
characteristics of the nonlinear dynamics. The flexibility effect can be verified by 
comparing the rigid body dynamics and the flexible body dynamics. The config- 
uration change effect and Coriolis and centrifugal force effect can be verified by 
changing the amplitude and the frequency of a sinusoidal excitation. In the simula- 
tions, a RALF model is tested with 30 LB payload in order to magnify the effects 
of the nonlinearity. In Fig. 5.4 and 5.5, several time responses of the rigid RALF 
model are compared by changing the amplitude (0.05 rad, 0.1 rad, 0.5 rad) with 
the fixed frequency (1.5 Hz) and the frequency (1.5 Hz, 4 Hz, 7 Hz) with the fixed 
amplitude (0.05 rad) of a sinusoidal excitation. Time responses are measured by 
several sensors - encoder, tachometer, and accelerometer. The TREETOPS model 
is used because the tip acceleration is easily measured using the given sensor option. 
The differences between each case are not clear in the time response. However, the 
phase plots in Fig. 5.6, 5.7 and 5.S are another method to detect the nonlinearity. 
The ellipse of harmonic excitation is distorted by the nonlinearity as shown in Fig. 
5.7 where the sinusoidal responses of different amplitude are drawn in the same 
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Fig. 5.4: Time response of the rigid RALF model 

for the different excitation amplitudes 
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Fig. 5.5: Time response of the rigid RALF model 

for the different excitation frequencies 
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plot. The acceleration at the tip shows the difference even in the time response 
as shown in Fig. 5.4. FFT of the tip acceleration is another method to check the 
nonlinearity. Fig. 5.9 shows the power spectral density of the tip acceleration in 
response to the different amplitudes and frequencies of the sinusoidal excitation. 
To obtain the sharp peak, 4096 data points between 0.9 second and 5.0 second are 
used. In order to reduce leakage, the Hanning window is used. The relative differ- 
ence between the first harmonic peak and the second harmonic peak decreases (a 
factor of more than SO) as the amplitude of sinusoidal excitation increases as shown 
in Table 5.1. a. However, the relative difference between the first harmonic peak and 
the second harmonic peak does not change significantly (a factor of less than 3) by 
changing the frequency of the sinusoidal excitation as shown in Table 5.1.b. 

Similar analysis has been done for the flexible body. Fig. 5.10 and Fig. 5.11 show 
the time responses of #i, and acceleration at the tip by changing the amplitude 
and the frequency of the sinusoidal excitation. In these figures, the distortion of the 
sinusoidal response can be detected even in the time response. Fig. 5.12 and Fig. 
5.13 show the phase plot of rotation angles and modal coordinates. Compared to 
the rigid body case, the effect of nonlinear dynamics is more significant as observed 
by changing either the amplitude or the frequency of the sinusoidal excitation. 
For instance, the circle of harmonic excitation with 0.5 rad amplitude and 7 Hz 
frequency is distorted as shown in Fig. 5.14. Fig. 5.15 shows the power spectral 
density of the acceleration at the tip. By increasing the amplitude of sinusoidal 
excitation, the broad range harmonic peaks including system natural frequencies 
are observable as shown in Fig. 5.15.a. By increasing the frequency of the sinusoidal 
excitation, the nonharmonic peaks occur between the harmonic peaks as shown in 
Fig. 5.15.b. Table 5.2.a and 5.2.b shows the first several harmonic peak values and 
their relative difference for different excitation amplitude and frequency respectively. 
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Fig. 5. 9. a: The power spectral density of the tip acceleration 

for the different excitation amplitudes 
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Fig. 5.9.b: The power spectral density of the tip acceleration 

for the different excitation frequencies 
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Table 5.1. a: 


0.05 rad 
0.1 rad 
0.5 rad 

Table 5.1.b: 


1.5 Hz 

4.0 Hz 

7.0 Hz 


Peaks of power spectrum density 

of rigid RALF model 

for different excitation amplitudes 
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Peaks of power spectrum density 
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Fig. 5.11: Time response of the flexible RALF model 

for the different excitation frequencies 
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Fig. 5. 15. a: The power spectral density of the tip acceleration 

for the different excitation amplitudes 



Frequency (Hz) 

Fig. 5.15.b: The power spectral density of the tip acceleration 

for the different excitation frequencies 


Table 5. 2. a: Peaks of power density spectrum 

of flexible RALF model 
for different excitation amplitudes 



1st 

2nd 

3rd 

lst/2nd 


harmonic 

harmonic 

hannonic 

harmonic 

0.05 rad 

1.73E5 

1.53E2 

3.36E2 

1.13E3 

0.1 rad 

7.00E5 

2.55E3 

1.22E3 

2.74E2 

0.5 rad 

2.47E7 

1.50E6 

2.96E5 

1.65E1 


Table 5.2.b: Peaks of power density spectrum 

flexible RALF model 
for different excitation frequencies 



1st 

2nd 

3rd 

1st /2nd 


harmonic 

harmonic 

harmonic 

harmonic 

1.5 Hz 

1.73E5 

1.53E2 

3.36E2 

1.13E3 

4.0 Hz 

7.25E6 

2.1SE5 

2.40E3 

3.33E1 

7.0 Hz 

1.SSE6 

2.3SE5 

1.76E4 

0.79E1 
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The significance of the harmonics grows dramatically in both cases (factor of about 
70 and 140) Therefore, the flexibility has a significant effect on the system dynamics 
as shown by comparison of Fig. 5.9 and Fig. 5.15. 

5.3.2 Experiment of Nonlinear Dynamics 

To test the nonlinear dynamics, fast motion of each link is required. The flow 
rate of the current hydraulic servovalve is 5 gpm. The effective piston area is 
different for each cylinder (See Appendix C.2). Therefore, the maximum linear 

speeds of the lower link actuator are 6.127 in/sec and 8.170 in/sec for extension and 

retraction, respectively [28]. The maximum linear speeds of the upper link actuator 
are 2.320 in/sec and 3.268 in/sec for extension and retraction, respectively. The 
transformation of the hydraulic cylinder’s displacement to the joint angle rotation 
results in 18.8 deg/sec, 24.4 deg/sec for the lower link and 7.34 deg/sec, 10.29 
deg/sec for the upper link at the nominal joint angles of the experiment. The 
maximum displacement of the cylinder is calculated as follows. 

9 = ASinut (5-1) 

9 = —Smut (5-2) 

u> 

Therefore, the maximum displacement of the lower link with 1.5 Hz frequency is 
0.65 in (1.99 deg or 0.035 rad) and 0.87 in (2.59 deg or 0.045 rad) for extension and 
retraction. The maximum displacement of the upper link with 1.5 Hz frequency is 
0.25 in (0.78 deg or 0.014 rad) and 0.35 in (1.09 deg or 0.019 rad) for extension 
and retraction. As we can see from comparison with simulation, the speed is too 
small for testing the nonlinear dynamics. Within the current maximum speed, the 
power spectral density is measured for different amplitudes and frequencies. Fig. 
5.16.a is the power spectral density by increasing the amplitude from 2 in, to 4 
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in, to 6.5 in of the lower link and 2.5 in of the upper link with 0.15 Hz frequency. 
The relative peaks decrease. Fig. 5.16.b is the power spectral density obtained by 
increasing the frequency from 0.36 Hz. to 0.59 Hz. to 0.£>2 Hz with 1 in amplitude of 
both actuators. Fig. 5.16.C is the power spectral density obtained by increasing the 
frequency from 3.6 Hz, to 5.5 Hz, to 7.9 Hz with 0.1 in amplitude of both actuators. 
As shown in these figures, the harmonic peaks grow slightly in the wide range with 
increasing excitation amplitude and frequency. A slight nonlinearity is observed 
under the current hydraulic cylinder's speed. 

5.4 Hydraulic Actuator Dynamics 

In this section, the hydraulic actuator dynamics are modelled and compared 
with dvnamics of the more commonly used electric D.C. motor. The difference in 
velocity feedback for the two actuators is shown to result in significantly different 
joint behavior for the two cases for a flexible structure. A simple single link system 
is used to illustrate this point. 

Dynamics of an asymmetric or a single - rod hydraulic cylinder is more compli- 
cated than that of a symmetric or a double - rod hydraulic cylinder. Because the 
area of each chamber of the cylinder is different, the retraction and extension speed 
is different.. Therefore, the dynamics of the asymmetric hydraulic cylinder is non- 
linear, As shown in many studies on an asymmetric hydraulic cylinder [24,75,36], 
the exact dynamics of a real actuator is high order and therefore hard to apply in 
real time control. Therefore, in this thesis, a linear time invariant model is used 
bv assuming the actuator is a symmetric hydraulic motor. Furthermore, in a crank 
mechanism, the relationship between the link rotation angle and the hydraulic pis- 
ton displacement is nonlinear so that the dynamics change depend on the operating 
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point. 

As shown in Appendix B, the governing equations of an electric armature con- 
trolled servomotor and a hydraulic motor are similar each other. In the electric 
armature controlled servomotor, the dynamic equations are 

T = Kii a = JO (5.3) 

I o % + Ral* = fia - Kj (5.4) 

at 

In the hydraulic motor, the dynamic equations are 

T = D m Pi = J6 (5.5) 

+ KctP = K q X v - D m 0 (5-6) 

4 O at 

Equations (5.3) and (5.4) axe similar to equations (5.5) and (5.6), respectively. The 
block diagram of the two models also show the similarity as shown in Fig. 5. IT. a 
and S.lT.b. The actuator dynamics can be embedded into the plant dynamics. The 
linearized state equations of the plant are 

ip = ApXp + BpUp (5.7) 

y = CpX p (5.8) 

The state equation of the actuator is 

x„ = -4 a x a + B a u a (5.9) 

Up = c a x a (5.10) 

If the loop is closed, the input vector is 

U a = G a Xr CpXp 


(5.11) 
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k e electric armature controlled servomotor 
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Fig. 5 . 17 .b: The block diagram of the hydraulic motor 
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Fig. 5.18 clarifies the interconnection of the plant and the controller. The corre- 
sponding augmented state equations axe 


\ 

ip 


-4 p BpC a 


fxj 

0 


> — 
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x a 
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-B a C p A a 


X a J 

B a Ga 


It 


y = CpXp 


(5.12) 

(5.13) 


This general form of augmented equation can apply to a flexible robot. The state 
equations of a flexible robot are 
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(5.14) 


where 


M = 


M T r M rS 
Mfr Mff 



(5.15) 


M rr , M/f, and M r f are the generalized inertia matrices which are related to the 
rigid body, the flexible body, and the coupling between the rigid body and the 
flexible body respectively. Kff is the generalized stiffness matrix. 9 and £ are 
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Fig. 5.18: The block diagram interconnection between the plant and the actuator 
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the generalized coordinates of the rigid body motion and the flexible bod\ motion 
respectively. 

Using the actuator equations, the state equations of the flexible robot augmented 
with the electric armature controlled servomotor and the hydraulic motor are as 
follows. 
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(5.16) 


(5.17) 


The only difference in the models is the values of each parameter as shown in 
Table 5.3 [3,59]. As shown in the block diagrams of Fig. 5.17, the velocity is fed 
back when the actuator dynamics is included. If the velocity feedback terms, ^ 
and £m - ar e n ot included, the system poles are determined by the eigenvalues of 
($/ _ A p )(sl — A c ) as shown in Eq. 5.12. The poles of the combined system are 
the combination of those of A v and A a . However, if the velocity feedback terms 
are included, the system poles are changed. Fig. 5.19 shows the root locus of a 
single flexible link when the velocity feedback terms of the electric motor and the 
hydraulic motor are changed from zero to the values shown in Table 5.3. When 
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Table 5.3: Comparison of system parameters between 

a electrical motor and a hydraulic motor 
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the electric motor is added, the closed loop system poles move a little bit from the 
open loop system poles. On the other hand, when the hydraulic motor is added, 
the closed loop system poles move almost to the open loop system zeros. This is 
due to the large value of compared to 

This phenomenon also can be observed in the Bode plot. Fig. 5.20 and 5.21 show 
the Bode plot for the angular position of the single flexible link with the electric 
motor and the hydraulic motor respectively. The several peaks observed in Fig. 5.20 
can not be found in Fig. 5.21 because the system pole is located near the system 
zero. The closer a pole is to a zero, the less that mode appears in the output. 
In the experiments, Fig. 5. 22. a is the Bode plot for the angular velocity of the 
single flexible link when it is excited by the electric motor. Because of the encoder 
resolution, the angular velocity is measured. Fig. 5.22.b shows the Bode plots for 
the positions of two hydraulic cylinders of RALF. The peaks of corresponding to the 
svstem frequencies are not found as predicted in Fig. 5.21. Note that the different 
size actuators are compared because of the hardware limitations. Therefore, only 
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5. 19. a: The root locus of a single flexible link with the electric motor 
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Fig. 5.19.b: The root locus of a single flexible link with the hydraulic motor 
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Fig. 5. 22. a: Measured Bode plot of a single flexible link 

with the electric motor for the angular velocity to input 
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Fig. 5.22.b: Measured Bode plot of each link of RALF 

with the hydraulic cylinder for the cylinder position to input 
(a) the lower link, (b) the upper link 
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the qualitative characteristic can be compared. From the above discussion, the 
actuator dynamics has a significant effect on the flexible dynamics. The pole-zero 
cancellation observed for joint variables as the output is not observed for other 
outputs, of course. Tip position or strain will continue to show the resonances of 
the arm. In addition, it should be noted that high velocity feedback gains in the 
control will also increase the coupling term and thus make the actuator stiff with 
respect to the disturbances caused by arm dynamics. 

5.5 Conclusion 

The TREETOPS model of RALF is used for the verification of the analytical 
model. The step responses of both models match each other very well. It is ■verified 
that SVD is a stable algorithm for a constrained dynamic system simulation. For 
verification of the nonlinear dynamics of RALF, a sinusoidal excitation method is 
used with various amplitudes and frequencies. The response is measured at the tip 
acceleration. As shown in the phase plots and the power spectra, the nonlinear 
dynamics of a flexible robot is significant compared to that of a rigid robot. The 
differences the amplitude of between harmonic peaks as shown in the power spectra 
can be used to check the degree of nonlinearity. However, because of the speed 
limitation of the current hydraulic cylinder, the nonlinear dynamics of RALF is not 
fully checked experimentally. 

The actuator dynamics effect on the flexible robot is also investigated. Even 
though the electric motor and the hydraulic motor have the same form of dynamic 
equations, the effect on the flexible robot can be different due to the difference of 
their system parameters. The differences are analyzed by the root locus and the 
Bode plot, both theoretically and experimentally. 
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CHAPTER VI 


Control of RALF 


6.1 Introduction 

The current RALF structure has some complicated components and geometries 
which axe haxd to describe analytically. As mentioned in the previous chapter, the 
actuator dynamics is an important component in the flexible arm dynamics. The 
dynamics of an asymmetric hydraulic cylinder is nonlinear. In a large and fast mo- 
tion, the nonlinear effect of the asymmetric hydraulic cylinder might be significant. 
The crank mechanism of the hydraulic cylinder and the offsets ( See Appendix Fig. 
C.l) make the symbolic modelling difficult. Therefore, before implementing a mod- 
ern control algorithm which requires an accurate dynamic model, a single input 
multi output controller is implemented using cylinder position and beam strain at 
each link independently. The performance of the classical controller can be used as 
the base performance for implementation of a more advanced control algorithm in 

the future. 

6.2 Controller Design and Experiment 

The position of each hydraulic cylinder is fed back using an LVDT sensor for the 
rigid body motion control. Because the hydraulic actuator has a velocity feedback 
loop internally, the position feedback with lag compensator has been generally used 
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for the position control of the hydraulic system [49]. 


u, 3 -r 

^ 1 > 7T 

s T j 


for i = 1,2 


(6.1) 


The lag compensator Eq. (6.1) is used to compensate the steady state error due 
to the friction. However, because the gravitational torque is much larger than the 
dynamic error correction torque, the compensation speed is slow. Therefore, the 
addtional constant torque is included to compensate the external torque due to 
gravity force. Gains and coefficients of the lag compensator are adjusted to obtain 
the underdamped step response. The parameters of the lag compensator are shown 
in Table 6.1. A cycloid curve is used for the trajectory profile. The reference tra- 
jectories used consist of the combination of the extension and retraction motion of 
two joints as shown in Table 6.2. The controller is digitized using a bilinear trans- 
formation method, s = ^7^7, and implemented in a Micro- Vax II computer. The 
sampling frequency used, T s , is S ms. The controller implementation is described 
in Appendix C. The control loop of the upper link is closed first and followed by 
the lower link. Fig. 6.1 and 6.2 are experimentally measured time histories of each 
cylinder position in cases of small motion and large motion respectively. In the two 
cases, the system follows the reference trajectories well. Fig. 6.3 shows measured 
time history of each cylinder position when a 30 lb payload is added at the tip. The 
system still follows the reference trajectories well. 

The strain is fed back to suppress the beam vibration. The strain feedback 
controller can be designed by root locus. The direct strain feedback pushes system 
poles into the right half plane as shown in Fig. 6.4. Therefore, the strain rate 
feedback is required to give flexible mode damping. The strain rate can be obtained 
by numerical differentiation of the measured strain. However, a low pass filter is 
additionally needed to reduce the noise effect. The strain rate also can be obtained 


Table 6.1: Parameters of a lag compensator 




0 
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Lower Link 0.5 Hz 

2.0 

5000. 


Upper Link 2.0 Hz 

4.0 

5000. 


Table 6. 2. a: Reference trajectory I 

(small 

motion) 


1st 

2nd 

3rd 

4th 

Lower Link [in]* 

2.0 

2.0 

4.0 

2.0 

Upper Link [in]* 

2.0 

4.0 

4.0 

* 

2.0 

Moving Time [sec] 

2.0 

2.0 

2.0 

2.0 

Pause Time [sec] 

10. 

10. 

10. 

10. 

Table 6.2.b: Reference trajectory II 

(large 

motion) 


1st 

2nd 

3rd 

4th 

Lower Link [in]" 

2.0 

2.0 

8.0 

2.0 

Upper Link [in]’ 

2.0 

4.0 

4.0 

2.0 

Moving Time [sec] 

2.0 

1.0 

1.0 

1.0 

Pause Time [sec] 

10. 

15. 

15. 

15. 


* absolute displacement of the cylinder 
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Fig. G.l.a: Measured 
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Fig. 6.2.b: Measured time response of RALF to a large step input 

a) Upper link strain, b) Lower link strain 
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Fig. 6.3.b: Measured time response of RALF to a large step input 

when a payload is attached 
a) Upper link strain, b) Lower link strain 
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from a filter. The generalized structural filter presented by Wie [78] is an extension 
of the phase -lead, -lag, bandpass, and notch filter. This approach for controlling 
flexible systems has been proposed for control of large space structures. We will 
explore its use for flexible arm control. The concept is based on various pole-zero 
patterns that can be realized from a second order filter presented as follows. 

(s 2 /^ + 2Cz^M + l) / 6 o) 

(sV^ 2 + 2Cp5/u;p + l) 

By different choices of the coefficients of the above second order filter, several fre- 
quency shaping filters such as phase lead and lag, notch, bandpass can be easily 
realized. In addition to these minimum- phase filters, various nonminimum- phase 
filters can also be realized from this second order filter. For stabilization of flexible 
modes, a phase stabilization technique and a gain stabilization technique using the 
phase lead or the notch filter have been used in practice [77,17,81]. Phase stabi- 
lization provides the proper gain and phase characteristics at the desired frequency 
to obtain a closed loop damping. Gain stabilization provides attenuation of the 
control loop gain at the desired frequency to ensure stability. In the system whose 
parameters are not known precisely, phase stabilization is more desirable. High loop 
gain can be realized if the desired mode does not cross the -180 deg. i n*360 deg. 

The first approach for vibration suppression is to change the phase by the phase 
lead filter. The phase lead filter can be obtained by setting > w z and C P = G = Cc 
as shown in Fig. 6.5.a. As shown in Fig. 6.5.b, the gain increase at higher frequency 
can be approximated as 40 * Ratio ^ greater than 2 should be avoided 

because a large ratio amplifies measurement noise. The filter coefficients are chosen 
by setting lj p and u ) z on either side of the mode to be stabilized. The damping ratio 
of the filter can be changed arbitrarily. For Cc = 1, the conventional phase lead filter 
with poles and zeros on the reed axis can be realized. Fig. 6.6 and 6.7 show the 
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root locus when the phase lead filter is used. The filter gain is varied. The phase 
lead filter can stabilize the modes of the upper link and the lower link conditionally. 
The effect of filter damping Q is not significant. 

The second approach is to employ the nonminimum phase allpass filter which 
maintains the control loop gain and provides the proper phase lag to the flexible 
mode as shown in Fig. 6.S.b. As shown in Fig. 6.S.a, the nonminimum phase 
allpass filter can be obtained by setting oj p = u: z = u; c and C P = “G* Fig. 6.9, 
several filter frequencies, ui p and are tried for the stabilization of the upper link 
mode. 30 filter gains which are logarithmically equally spaced between 1 and 1000 
are chosen. If the filter frequency is higher than that of the upper link mode to 
be stabilized as shown in Fig. 6.9.C, it is hard to stabilize the upper link mode by 
the filter. However, if filter frequency is lower than that of the upper link mode, 
the filter can stabilize the upper link mode. As shown in Fig. 6. 9. a and 6.9.b, 
the system poles move more slowly by the lower filter frequency for the same gain. 
Therefore, the lower filter frequency of Fig. 6. 9. a is more desirable. In Fig. 6.10, 
the same technique is applied to the lower link mode. In this case, if the filter 
frequency is lower than that of the lower link mode to be stabilized, it is hard to 
stabilize the lower link mode by the filter as shown in Fig. 6. 10. a. However, if the 
filter frequency is higher than that of the lower link mode, -the filter can stabilize the 
lower link mode conditionally. As shown in Fig. 6.10.b and 6.10.c, the system poles 
also move more slowly by the lower filter frequency for the same gain. Therefore, 
the filter whose frequency is close to the lower link mode is more desirable. Fig. 6.11 
and 6.12 are simulated time responses of an analytical model to verify the results 
observed in Fig. 6.9 and 6.10 respectively. As shown in Fig. 6.9 and Fig. 6.10, the 
system response is dominated by the lower link mode because the lower link mode 
is close to the imaginary axis. As expected in the previous discussion. Fig. 6, 11. a 
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Fig. 6.5.a: Pole - Zero plot of the phase lead filter 




Fig. 6.5.b: Gain - Phase plot of the phase lead filter 
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Fig. 6.7: Root locus of the lower link when strain is fed back 

through a phase lead filter 
a) Cc = lv b)Cc=-7 





Fig. 6.8.a: Pole - Zero plot of nonminimum-phase allpass filter 



Fig. 6.8.b: Gain - Phase plot of nonminimum-phase allpass filter 




Fig. 6.9: The root locus of the upper link when the strain is fed back 

through a nonminimum phase allpass filter 
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Fig. 6.10: The root locus of the lower link when the strain is fed back 

through a nonminimum phase allpass filter 
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and Fig. 6.12.b show the best response for the same filter gain (K f - 10). 

In experiment, the strain is detected near the connecting points of the hydraulic 
cylinder and the link (Appendix Fig. 6.11). Fig. 6.1.b, 6.2.b. and 6.3.b show 
the strain histories of the previous three cases. The beam does not vibrate much 
except when starting the motion because of current speed limit and smooth cycloid 
trajectory. In order to magnify the vibration, a faster reference trajectory is chosen. 
New reference trajectory is 2 in move of the lower beam and 1 in move of the upper 
link during 0.1 sec. The controller is also digitized using a bilinear transformation 

method. 

First, the strain is fed back through the phase lead filter. Fig. 6.13.a shows the 
strain history of the upper link strain when no strain is fed back. As shown in Fig. 
6.13.b, the frequency spectrum has two peaks corresponding to the lower and the 
upper link frequency respectively. Fig. 6. 14. a shows the strain history of the upper 
link when the strain is fed back by the phase lead filter. As shown in Fig. 6.14.b, 
the peak corresponding to the upper link is significantly reduced. Fig. 6.15.a shows 
the strain history of the lower link strain when no strain is fed back. As shown in 
Fig. 6.15.b, frequency spectrum has two peaks corresponding to the lower and the 
upper link frequency respectively. Fig. 6.16.a shows the strain history of the lower 
link when the strain is fed back by the phase lead filter. As shown m Fig. 6.16.b, 
the peak corresponding to the lower link is also reduced. However, a higher gam 
makes the system unstable. Because the phase lead filter excites higher modes while 
changing the phase, high gain may excite the unmodelled dynamics. Furthermore, 
the restriction on & limits the selection of filter parameters. 

Secondly, the strain is fed back through the nonminimum-phase allpass filter to 
give a proper phase margin without changing the gain. Because this filter does not 
change the loop gain, it is easy to adjust the filter coefficients and high feedback 
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gain is possible. Based on the root locus analysis, the filter frequency for the upper 
link is chosen as 3.0 Hz which is below the upper link frequency (9.S Hz). The filter 
frequency for the lower link is chosen as 5.4 Hz which is the same as the lower link 
frequency (5.4 Hz). The filter damping ratio is chosen as 0.7. Therefore, coefficients 
of the nonminimum phase allpass filter used are shown in Table 6.3. 

Fig. 6.17 and 6.18 show the strain histories and frequency spectrum of the 
upper link and the lower link when the strain is fed back through the nonminimum 
phase allpass filter. The filter is designed for each link using the information of 
the corresponding strain. The vibration is significantly reduced. However, the 
frequency corresponding to the other link still remains in the time response as 
shown in Fig. 6. 17. a and 6. 18. a. Suppression of both frequencies at the same time 
is difficult because it requires a higher order filter. Furthermore, because the two 
frequencies have different phase, the phase adjustment for one frequency makes the 
other frequency unstable. 

Until now, each controller is designed independently when only the correspond- 
ing link is moved. Next, the responses are measured when the both joints axe moved 
at the same time. Fig. 6. 19. a shows the time response of the upper link when no 
strain is fed back. Fig. 6.19.b shows the corresponding frequency spectrum. How- 
ever, when the strain is fed back through the nonminimum phase allpass filter, the 
vibration is significantly reduced as shown in Fig. 6.20. Similarly, the vibration of 
the lower link also reduced as shown in Fig. 6.21 and Fig. 6.22. Compared to Fig. 
6. 11. a and Fig. 6.12.b, time responses of each link have more damping because of 
the structural damping. Furthermore, measured strain has offset due to the gravity 
effect . 

In order to check the robustness to frequency change, the 30 lb payload is at- 
tached at the tip. The same type test as above has been performed. Fig. 6.23 and 
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6.24 show the strain histories of the upper and the lower link respectively when no 
strain is fed back. By attaching the payload, the frequency of the upper link be- 
comes the fundamental frequency. That frequency is dominant as shown in the both 
frequency spectra. Fig. 6.25 and 6.26 show the strain histories of the two links when 
the strain is fed back through the nonminimum phase allpass filter. The original 
filter makes the system unstable. Even though the gain is reduced to one fourth of 
filter gain, the strain feedback cannot reduce the vibration. However, by changing 
the filter frequencies of the upper link slightly, the vibrations are reduced as shown 
in Fig. 6.27 and 6.28. The filter coefficients of the upper link are changed as shown 
in Table 6.4. However, the better response is hard to obtain just by the filtering 
technique because the two system frequencies are close together when the payload 
is attached. Fig. 6.29, 6.30. and 6.31 show the strain histories when strain is fed 
back through the nonminimum phase allpass filter. These use the same trajectory 
as Fig. 6.1. 6.2, and 6.3. 


6.3 Conclusion 

The independent control scheme using position and strain is applied for the 
control of RALF. Position feedback with a lag compensator is successfully imple- 
mented for the rigid body motion control. However, the direct strain feedback does 
not reduce the vibration. For the vibration suppression, a phase lead and a non- 
minimum phase allpass filter are used. In experiment, a nonminimum phase allpass 
filter is easier to implement than a phase lead filter. As shown in strain histories 
of an analytical model and experiment, the proper filter frequencies can be selected 
using root locus. Therefore, even though there is a limitation in using these filters 
for a system with changing frequencies, the beam vibration is reduced significantly 
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Table 6.3: Nonminimum phase allpass filter coefficients 



Lower 

Upper 

Gain 

50. 

SO. 


5.4 

3. 

uj z 

5.4 

3. 

Cp 

0.7 

0.7 

G 

-0.7 

-0.7 


Table 6.4: Nonminimum phase allpass filter coefficients 

when the payload is attached 



Lower 

Upper 

Gain 

12. 

20. 


5.4 

1. 


5.4 

1. 

Cp 

0.7 

0.7 

G 

-0.7 

-0.7 


by the nonminimum phase allpass filter. The sensitivity to payload variation is 
pronounced, however and is a deterrent to using the allpass filter in manipulator 


control. 
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Fig. 6.14: Measured time history and frequency spectrum 

of the upper link strain 

when strain is fed back through a phase-lead niter 
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Fig. 6.16: Measured time history and frequency spectrum 

of the lower link strain 

when strain is fed back through a phase-lead filter 
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Fig. 6.19: Measured time history and frequency spectrum 

of the upper link strain 
when no strain is fed back 
(both joints move) 
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Fig. 6.20: Measured time history and frequency spectrum 

of the upper link strain 

when strain is fed back through a no nmin i m um- phase allpass niter 
(both joints move) 
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Fig. 6.23: Measured time history and frequency spectrum 

of the upper link strain 
when no strain is fed back 
(both joints move & payload) 
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Fig. 6.24: Measured time history and frequen 

of the upper link strain 
when strain is fed back through a i 
(both joints move &z payload) 
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Fig. 6.25: Measured time history and frequency spectrum 

of the lower link strain 
when no strain is fed back 
(both joints move & payload) 
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Fig. 6.27: Measured time history and frequency spectrum 

of the upper link strain 

when strain is fed back through an adjusted nonminimum-phase 
allpass filter (both joints move & payload) 
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Fig. 6.28: Measured time history and frequency spectrum 

of the lower link strain 

when strain is fed back through an adjusted nonminimum-phase 
allpass filter (both joints move &z payload) 
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Fig. 6.29: Measured time histories 

of the upper and the lower link strain 

when strain is fed back through a nonminimum- phase allpass filter 
(both joints move k small motion) 
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Fig- 6.30: Measured time histories 

of the upper and the lower link strain 

when strain is fed back through a nonminimum- phase allpass filter 
(both joints move L. large motion) 
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CHAPTER VII 


Summary and Future Work 


7.1 Summary 

For the purpose of exploring possible industrial applications, a large lightweight 
manipulator was previously constructed at the Flexible Automation Laboratory 
of the Georgia Institute of Technology. For added payload capacity, the flexible 
parallel mechanism is added to sustain the load without increasing weight very 
much. If RALF is assumed as rigid, the dynamics of RALF can be easily analyzed 
using the theories which have been developed for rigid robots. However, in a flexible 
robot, there axe several problems which are not observed in a rigid robot. In this 
thesis, solutions of these problems are evaluated. This section briefly summarizes 
the important observations and results of the previous sections. 

First, the dynamic equations of motion of flexible robots are complicated because 
of the link deformation which is expressed by the summation of admissible functions. 
The existing symbolic method generates complicated equations of motion. This 
thesis develops a structurally well organized and computationally efficient form of 
the equations of motion using the Jacobian matrix. The Jacobian matrices are 
derived from the position vectors directly using the MJac function of SMP. Using 
this Jacobian matrix, the mass matrices and the gravity force vectors are derived. 
The velocity coupling terms are derived from the mass matrices using the Christoffel 
symbol. One problem of this method is that the velocity coupling terms comprise 
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many elements. However, using some structural properties, such as symmetry and 
reflective coupling, these elements can be expressed by a few parameters. The 
resulting equations of motion have a structure which is useful to reduce the number 
of terms to be calculated, to check correctness, or to extend the model to higher 

order. 

Secondly, the parallel link mechanism with flexible links requires nonlinear al- 
gebraic constraint equations to describe the relationship between angles within the 
closed loop chain. Therefore, a mixed set of differential equations and algebraic 
equations must be solved simultaneously for the constrained system dynamics. Be- 
cause the integration of these equations is sensitive to numerical error, the conven- 
tional numerical methods are not adequate to solve these equations. Among several 
methods available to solve these numerical difficulties, a coordinate partitioning 
method is used in this thesis. The generalized coordinates are transformed to a 
set of independent coordinates on the tangent hyperplane of the constraint surface. 
The transformation matrix used is the null space matrix of the constraint Jacobian 
matrix which can be derived by the singular value decomposition. Because the La- 
grange multiplier is eliminated from the equations of motion by the transformation, 
the constraint forces do not affect the motion on the constraint surface and only 
maintain the system on the constraint surface. Therefore, numerical integration 
on the tangent hyperplane does not result in significant constraint violation. This 
thesis applies this constrained system analysis technique to a flexible parallel link 
mechanism. It is numerically verified that the SVD method is a stable algorithm 
to solve a constrained system. 

Third, mode shape functions are employed to describe the elastic deformation of 
the flexible manipulator. In the assumed mode method, mode shape functions need 
only to be admissible functions which satisfy the geometric boundary conditions. 
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However, a large number of modes are required to obtain accurate frequencies. The 
number of modes to be included can be reduced by choosing appropriate functions 
which satisfy static equilibrium at the interface between the links. These functions 
can be derived systematically using component mode synthesis. Various component 
mode synthesis approaches have been developed depending on the assumed bound- 
ary condition. Different boundary condition assumptions require different coordi- 
nate systems to describe the elastic deformation. Therefore, the method which fits 
with the current coordinate system was chosen. In this thesis, proper mode shape 
functions for a reduced order model of RALF are obtained using the loaded interface 
component mode synthesis method. Comparison between the finite element model, 
the component mode synthesis model, and the analytical model of RALF show that 
the mode shape functions which axe determined by component mode synthesis im- 
prove the convergence. However, there are some discrepancies between the analytic 
model and the experiment because the real system has a complex structure which 
is difficult to analyze with the analytical method. For example, the upper link 
requires a detailed model to predict the better approximate boundary condition. 
These discrepancies axe explained using a simplified and a detailed finite element 
model in this thesis. 

Fourth, a direct comparison between the analytical model and the real system 
was difficult because of the unmodelled dynamics and the parameter uncertainty of 
the real system. Therefore, as an intermediate step, a general multibody flexible 
dynamics code - TREETOPS - was used for the verification. The step response of 
the models match each other very well. 

Fifth, it is an important job for control to detect the degree of the nonlinearity in 
the dynamics. In this thesis, the nonlinearity of RALF is studied using a sinusoidal 
excitation. As shown in the phase plot and power spectra, the flexibility effect in 
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nonlinear dynamics is significant. However, the nonlinearity of RALF is not fully 
excited due to the speed limitation of the hydraulic cylinder. 

Sixth, the actuator dynamic effect on the flexible robot was investigated. Even 
though the electric motor and the hydraulic motor have the same form of dynamic 
equations, their effect on the flexible robot is different due to the difference of 
their system parameters. In this thesis, the effects of an actuator dynamics on 
a flexible robot are analyzed by root loci and Bode plots, both theoretically and 
experimentally. 

Seventh, in order to show the base performance for an advanced control, this 
thesis uses a decentralized control scheme using cylinder position and beam strain 
for trajectory following and vibration suppression. The rigid body motion is con- 
trolled by the position feedback using a lag compensator. The beam vibration is 
reduced significantly by strain feedback using the nonminimum phase allpass filter. 
Using these simple controllers, a good result is obtained. 

7.2 Future Work 

First, the proposed derivation method of equations of motion using Jacobian 
matrix has been applied only to a two joint planar system, RALF. Furthermore, 
some procedures for deriving the velocity coupling terms are not computerized. In 
the future, an even more computerized derivation method is required for multi-link 
flexible body dynamics. 

Secondly, one problem of a coordinate partitioning method is the preservation 
of continuity in the basis of the nullspace. Since the constraint Jacobian matrix 
is time-varying, the basis of the nullspace which forms a basis orthogonal to the 
tangent plane of the constraint surface is also time-varying. However, performing 
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singular value decomposition at each integration step would result m a tremendous 
increase in computing time and requires numerical integration algorithms which use 
information only from the current step. To overcome this problem, Liang and Lance 
[39] proposed a differentiable null space method. This method generates a set of 
independent coordinates that are on the constraint surface rather than only on the 
tangent hyperplane. Another solution is the modification of the constraint violation 
method. Chang and Nikravesh [16] applied an adaptive algorithm to stabilize the 
violation. Park and Haug [57] developed a hybrid numerical method which combines 
a constraint stabilization method and a generalized coordinate partitioning method. 
As we can see here, the subject of numerical methods for solving a mixed set of 
algebraic and differential equations is still an open research area. 

Third, mode shape functions of RALF were derived by hand calculation using 
component mode synthesis method from the output of MSC/PAL. An computer- 
ized program is required for different configurations of the structure. Furthermore, 
natural frequencies and modes of the system can be varied due to several effects - 
payloads, contact with the environment, joint position and velocity feedback gains 
[11]. The study of these uncertainties will require much work in the future. 

Fourth, nonlinear dynamics verification is a difficult job and has not been fully 
studied. Chaotic vibration theory might be useful to study the nonlinear vibration 

[50]. 

Fifth, as mentioned before, there are some uncertainties in the dynamic model 
of the flexible robot under the current theories. Furthermore, the system dynamics 
has a strong interaction with the controller output because it changes the boundary 
condition of the beam. Therefore, control schemes which are sensitive to the model 
uncertainty are difficult to apply. Furthermore, even the characteristics of the linear 
control of the multi-link flexible robot have not fully been studied. Root loci help us 
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to understand the control characteristics of the flexible robot [7]. The combination 
of state space and frequency domain techniques through a model update procedure 
has been shown effective in determining the feedback gams of the flexible robot [9]. 
Recently, a frequency domain analysis of a multi input and multi output system has 
been studied [44,20] in many areas. Especially, the Multivariable Frequency Domain 
Toolbox of Matlab [1 2] implements several frequency domain design techniques such 
as the Nyquist Array method, the Characteristic Locus method, the Quasi- Classical 
Design method, and the Multivariable Root Locus method. Using these computer 
aided design techniques, better studies on the linear control of the multi -link flexible 

robot can be expected. 
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appendix a 

Derivation of Equations of Motion using Jacobian Matrix 


In this appendix, a Lagrangian method is used to derive the equations of motion 
for a flexible robot. The Jaeobian matrices are used to derive the mass matrices and 
gravity force vectors. The coefficients of centrifugal and Coriolis force are denved 

from the mass matrices. 

The total kinetic energy of an elastic link can be written as 


T = 



rlr p p p A v dxj, 


(A.l) 


where b is the number of links, r, is the velocity vector of any point on the elastic link 
p, and p„ .4,, l„ are the density, the area, and the length of link p respectively. The 
velocity vector can be expressed by the Jacobian matrix and generalised velocity 

vectors q p . 


r p — JpQp 


(A.2) 


Substitute (A.2) into (A.l), 

T = 1 £ / (JpQpViJp&PpApdxp 

2 p=X Jo 

= \ E ( t J v J vPv A v dx v)iv 

~ P=1 

Equation (A. 3) can be written in a scalar form as 

- i=i j=l p - 1 


(A.3) 


(A.4) 
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where n is the number of the generalized coordinates, D ijp is ij element of the mass 
matrices JJJ pPp A p dx p for the pth link. If M, } is defined as 

p=l 

Equation (A. 4) can be written as follows 

T= 1 -±±M„m (AS) 

- .=1 J=1 


The potential energy due to gravity is 


U. 


b r l 

= 5/o 


9 Tr rPp A P dx ? 

p=l ' 


(A.7) 


where g is the 3 x 1 gravity acceleration vector. The potential energy due to elastic 
deformation is 

o.-\i < A - 3 > 

~ p=l v/u P 

where E is Young’s modulus of elasticity, and I is the area moment of inertia, u is 
the elastic deflection which can be expressed by m modes and modal coordinates. 


Up(x,f) = Y^ll>pk{x)Zpk(t) 
Jt=l 

Therefore, the elastic energy can be rewritten as (A.ll) 


(A. 9) 




*i 6 

(A. 10) 

- P = 1 /t=l 


where 



(A.ll) 

If the generalized coordinate is expressed as 


{ Sj ) = {«„u} r . j = 1 n 

(A. 12) 
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where n (= b + m) is the number of the generalized coordinates, 9 P is the rigid 
body motion coordinates of link p. £ pk is the k th modal coordinates of link p, the 
elastic energy can be rewritten by the generalized coordinates as follows. 




1=1 j= i 


(A. 13) 


Using the kinetic energy and the potential energy, each term of Lagrange’s equa- 


tion 

d t dT dT d(Lg + U e ) _ ^ (A. 14) 

dt dq t dq, dq, 

can be written as follows. If the kinetic energy (A.6) is substituted into (A. 13), 


= -jr (E - E MijUj + E fa 

dvdqi dt " J=l J=i ai 


dMi j . 




(A. 15) 


where 


the first term becomes 


dMij _ A dM^dqk _ A Wj, 
dt ~ dq k dt dq k 9k 


dM i: , . . 


(A. 16) 


*S> = 


and the second term becomes 


I 1 = 

U<lt - ; = 1 Jt=l 

A » 1 dM, k . . 

= SSs-ST®* 


(A. IT) 


Next, if the potential energy (A.7) and (A.ll) are substituted into (A.13), the third 
term becomes 

dU, 


BU, _ V t" a T ^o ‘ 

= £/ o 'V 4 VV* 


Apdip 


(A. 18) 



where is the i th column of Jacobian matrix J p . The third term can be written 
by a scalar because the gravity is acting in one direction. Therefore, if the gravity 
vector has nonzero term in r th row, 

U g = G ig ( A - 19 ) 


where 


6 f ip 

Gj = > , / Jp{ r i i]Pp-^-pd x t 

p=i Jo 


(A. 20) 


where 7 p [r, z] is r th row and i th column of J p . The fourth term becomes 


dUe 

dq : 


1 n m 
U( b - i=i j=i 


= £ K 'in 


(A. 21) 


J=1 


Therefore, the Lagrangian equations of motion can be written symbolically as fol- 
lows. 


Y + £ Kaqj+ 

;=i j=i 

^ ^ 1 „ 5M ti , 5M.it dM }k ^ . . n _ 
£ £ ^(-s- 2 + — —)W + G.flr = rt 


(A. 22) 


U ti 2 % dt H 


dqi 


or 

^ Mij?j + 52 Kirti + £ £ Gjk(i)qjqk + G.s = t, (A.23) 

j=i j=l i=i *=i 

where q is the vector of generalized coordinates, M is the generalized mass matrix, 
K is the elastic stiffness matrix, C is the coefficient matrix of Coriolis and centrifugal 
forces, G is the gravity force vector, r is the generalized force vector. 
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APPENDIX B 

ACTUATOR DYNAMICS 


The dynamics of a hydraulic motor is similar to that of an armature-controlled 
servomotor. In this appendix, the dynamics of the electric motor and the hydraulic 
motor are summarized and compared. 

B.l Armature Controlled Servomotor 

Consider the armature-controlled dc motor shown in Fig. B.l [55]. In this 
system, For a constant field current, the flux becomes constant, and the torque 
becomes directly proportional to the armature current. 

T = Kii a ( B - 1 ) 

where K { is a motor- torque constant. For a constant flux, the induced voltage e b is 
directly proportional to the angular velocity 6. 

«-Kj> ( B -2) 

where K b is a back emf constant. The differential equation for the armature circuit 
is 

L a ^+R a i a + e b = e a (B.3) 

ut 

L , ^ + R.i. = Kj 

dt 


or 


(B.4) 



Ra — armature-winding resistance, ohms 

L a = armature-winding inductance, henrys 

i a = armature-winding current, amperes 

ij = field current, amperes 

e a = applied armature voltage, volts 

ej, = back emf, volts 

9 = angle of the motor, radians 

T = torque delivered by the motor, N - m 

J as equivalent moment of inertia of the motor and load 
referred to the motor shaft, Kg — m 2 

The armature current produces the torque which is applied to the inertia. 

T = Kii a = JO (B.5) 


B.2 Hydraulic Actuator 

Consider the servovalve controlled hydraulic motor shown in Fig. B.2 [49]. In 
this system, The linearized servovalve flow equations are 

Ql = K q x v - K c P l (B.6) 

where 

Q l = = load flow, m 3 /sec 

P L = P x — P 2 = load pressure difference, N/m 2 

From the continuity equation of each motor chamber, the continuity equation 


for all hydraulic actuators is 
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Q i, Q-i 


PuPt 


K q 

I<c 


V t 


0 

T 

J 


forward and return flows, m 3 /sec 

forward and return pressures, N/m 2 

valve displacement from neutral, m 

valve flow gain, m 3 / sec/m 

valve flow-pressure coefficient, m 4 . sec/ Kg 

volumetric displacement of motor, rr? /rad 

total leakage coefficient of motor, m 4 . sec/ Kg 

total contained volume of both chamber, m 3 

effective bulk modulus, N/m 2 

torque generated by motor, N — m 

total moment of inertia of motor and load 

referred to motor shaft, Kg — m} 


K ce = K c + Ctm 

= total flow-pressure coefficient, m 4 . sec/ Kg 

Thus the load flow Q L is consumed by flow to displace the actuator, leakage, and 
flow stored due to compressibility. The leakage is usually neglected. From (B.6) 

and (B.7), _ 

Dj + = K,x. - K„Pl (B.S) 

4 0 

or 

^Pl + K C 'Pl = K q x v - D m 6 (B .9) 


where 

The torque balance equation is 


T = P L D m = JO 


(B.10) 
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Equations (B.4) and (B.o) resemble equations (B.9) and (B.10), respectively. 
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APPENDIX C 


Controller Implementation 


C.l Software Implementation 

This section gives some comment on data acquisition, transducer voltage conver- 
sion, and the computer implementation of path planning of the current controller. 
The detailed descriptions of controller implementation were explained in J.D. Hug- 
gins’s thesis [28]. 

1) A set of assembly language subroutines that could be called from Fortran 
or C programs axe used to derive a programmable clock and the A/D and D/A 
board. The data is read into the buffer using the subroutine DTSBR and DTSBWB. 
DTSBR reads the A/D channels. DTSBWB causes the computer to wait until the 
buffer is full before processing. After initialization, the DTSBR and DTSBWB 
subroutines can be called repeatedly to transfer the data. 

2) Calculate the length of the actuators from read data. The minimum and 
maximum values of each cylinder length and its digital numbers are measured as 
in Table C.l. Therefore, the relation between the cylinder length y, and its digital 
numbers x, can be described as equation (C.l). 


Vi — m;x, + k (i = 1,2) 


(C.l) 



1S1 


The coefficient m, and 6 , can be calculated as follows 


mi 


28.750 - 34.406 
1080 + 1941 


= -1.4144 x 10 -3 


bi = 28.750 + 1.41444 x 10 -3 x 1080 = 31.566 


(C. 2 ) 

(C.3) 


m 2 


29.625 — 35.375 
1080 + 1941 


= -1.903343 x 10 ~ 3 


b 2 = 29.625 + 1.903343 x 10 ' 3 x 1080 = 31.6806 


(C.4) 

(C.5) 


3) As shown in Fig. C.l, the attachment points of actuators axe offset from the 
centerline of the beams. The initial angles of each link are calculated from geometry 
as shown in Fig. C. 2 . 


$io = $n — $12 + $13 — 36.5(5° (C.6) 

$20 = $21 — $22 + $23 “ (0.334 (C. 1 ) 

where the connecting link and the upper link are not parallel each other due to the 
bracket offset of the upper link. Therefore, $03 is the angle difference between the 
connecting link and the upper link. 

Table C.l: Cylinder length and its digital number 


Cylinder 


Length 

Number 

1 

Min 

28.750 in 

+ 1991 


Max 

34.406 in 

- 2008 

2 

Min 

29.625 in 

+ 10 S 0 


Max 35.375 in - 1941 
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4) The path planning algorithm needs the initial position, final position, and 
the time for the move. The flow chart of path planning and controller is depicted 

in Fig. C.3. 

5) The cylinder position is inversely proportional to the cylinder bore size. Be- 
cause the two cylinder bore diameters are different, the gains of two amplifiers are 
different from each other. 

C.2 Equipment List 
COMPUTER 

Model : Microvax II (vector of VS21W - A2) 

Company : Digital Equipment Co. 

ADDITIONAL BOARDS for Micro Vax 

Model : DT2769 - Real Time Clock Board 

Model : DT2785 - Analog I/O system 

8 channels A/D multiplexed 
2 channels D/A 
12 bit resolution 

Company '■ Data Translation 

SIGNAL ANALYZER and DISK STORAGE UNIT 
Model : 3562A - Digital Signal Analyzer 

Model : 9122 - Disk Storage Unit 

Company : Hewlett - Packard Co. 



1S3 


YO 



Figure C.l: Attachment points of actuators and bracket offsets of links 

















STRAIN GAUGE 

Model : A - 13 -250MQ - 350 

Company : Measurements Group Inc, Micro Measurement Div. 

STRAIN GAUGE AMPLIFIER 

Model : 3B18 - Wide band 

Company : Analog Devices 

HYDRAULIC COMPONENTS 

Servo Valves 

Model : 73 - 102A Two Stage Servovalves - 5 gpm 

Company : Moog, Inc 

Cylinder of Lower Link 

Model : H - PB - 2 Cylinder 

Bore : 2.0 in. 

Stroke : 20 in. 

Rod Diameter : 1.00 in. 

Seals : Teflon 

Weight • 35 lbs. 

Company : Atlas Cylinder Corp. 

Cylinders of Upper Link 

Model : N2C - 3.25 x 40 Cylinder 


Bore 


: 3.25 in. 

Stroke : 40 in. (modified to 17 in.) 

Rod Diameter : 1.75 in. 

Seals : Buna - N 

Weight : 52 lbs. 

Company : Hydroline Mfg. Co. 

POWER AMPLIFIER 

Model : BOP15 - 20M 

Spec. : 20amp, 15 volt, 300kHz crossover 

Company : KEPCO 
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APPENDIX D 
SMP Code 

In this appendix. SMP codes for the derivation of symbolic equations of motion 
of the lower and the upper links are described. Figure D.l is the flow chart of SMP 
code. In these codes, program is separated into several files because one program 
depletes computer memory. First, Jacobian matrices are derived from the position 
vectors using MJac function of SMP as shown in Fig. D.2. Several simplication 
procedures axe required to generate a compact form of output. Secondly, elements 
of the mass matrices and the gravity force vectors are obtained using Jacobian 
matrices as shown in Fig. D.3. Third, the symmetric elements and the reflective 
coupling elements of Christoffel symbol are generated as shown in Fig. D.4 and 
D.5, respectively. Properties which are unique in flexible body dynamics are also 
included in Fig. D.4. For obtaining more reduced elements of velocity coupling 
matrices, additional simplications of the output are required by the supervision of 
the analyst. Finally, velocity coupling matrices are derived as shown in Fig. D.6. 
using the Christoffel symbol and the matrices are simplified from Fig. D.4 an D.5. 
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Figure D.l: Flow chart of symbolic program 
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/* symbolic Program of Lower and Upper Link */ 

cl: Cos [thl] 
sl:Sin[thl] 
cl2 : Cos [ thl+th2 ] 
s 12 : S in [ thl+ th2 ] 

/* Position vectors */ 

/* 2 modes per each link */ 

rl:{xl cl - pll qll si - pl2 ql2 sl\ 

,xl si + pll qll cl + pl2 ql2 cl) 
r2 : { 11 cl-plle qll sl-pl2e ql2 si \ 

+ X2 Cl2-p21 q21 Sl2-p22 q22 sl2 \ 

,11 sl+plle qll cl+pl2e ql2 cl \ 

+ x2 Sl2+p21 q21 Cl2+p22 q22 cl2) 

/* generalized coordinates */ 

q: { thl, th2 , qll, ql2 , q21, q22 } 

/* Substitutions */ 

sbO:Sin[$x] A 2->l-Cos[$x] A 2 
sbl:cl cl2->c2-sl sl2 
sb2:cl sl2->s2+cl2 si 

sb5:pll qll + pl2 ql2 ->ul 
sb6:plle qll + pl2e ql2 ->ule 
sb7:p21 q21 + p22 q22 ->u2 

/* Find Jacobian matrices */ 

<XMJac 

jcl: Cb[MJac[rl, q] , {-cl, -si) ] /* Cb: Combine */ 

j c2 : Cb [MJac [ r2 , q] , {-cl,-cl2, -sl,-sl2} ] 

jvl:S[ jcl, { sb5 , sb6 , sb7 } ] /* S: Substitute */ 

j v2 : S [ j c2 , { sb5 , sb6 , sb7 } ] 

j vlt: Trams [ jvl] /* Trams: Tramspose */ 

jv2t : Trans [jv2] 


Figure D.2: Symbolic program of Jacobian matrices 
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/* Find Mass Matrices */ 


m: (al jvlt.jvl +■ a2 jv2t.jv2 ) 
mll:Ex[S[Ex[m[l,l]] , {sb0,sbl,sb2}] ] 

ml2:Ex[S[Ex[m[l,2]], {sb0,sbl,sb2}]] 

ml3 :Ex[S [Ex[m[l, 3 ] ] , { sbO , sbl, sb2 } ] J 
ml4:Ex[S [Ex[m[l, 4] ] , { sbO , sbl, sb2 } ] ] 
ml5:Ex[S[Ex[m[l,5] ] , ( sbO , sbl, sb2 } ] ] 
ml 6 : Ex [ S [ Ex [ m [ 1 , 6 ] ] , {sbO,sbl,sb2} ] ] 


/* Ex : Expension */ 


m22:Ex[S[Ex[in[2,2] ] , {sbO} ] ] 
m23 : Ex[S [Ex[m[2 , 3] ] , {sbO, sbl, sb2 } ] ] 
m24 :Ex[S [Ex[m[2 , 4] ] , {sb0,sbl,sb2} ] ] 
m25 :Ex[S [Ex[m[2 , 5] ] , {sbO } ] ] 
m26:Ex[S[Ex[m[2,6] ] , {sbO} ] ] 

m3 3 :Ex[S [Ex[m[3 , 3] ] , {sbO} ] ] 
m34 : Ex[S [Ex[m[3 , 4] ] , {sbO}]] 
m35 :Ex[S [ExCm[3 , 5] ] , {sb°,sbl,sb2} ] ] 
m36:Ex[S[Ex[m[3,6]] , {sb0,sbl,sb2} ] J 

m44:Ex[S[Ex[m[4,4] ] , {sbO}]] 
m45:Ex[S[Ex[m[4,5]] , ] 

m46 : Ex[S [Ex[m[4 , 6] ] , {sbO , sbl, sb2 } ] j 


m55 : Ex[S [Ex[m[5 , 5] ] , {sbO}]] 
m56:Ex[S[Ex[m[5, 6] ] , (sbO) ] ] 

m66:Ex[S[Ex[m[6,6] ] , {sbO}] ] 

/* Find Gravity Vectors */ 

glial jvl[2,l] + a2 jv2[2,l] 
g2:a2 jv2[2,2] 


Figure D.3: Symbolic program of mass matrices and gravity vectors 
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/* symmetric Property of Christoffel Symbol */ 

♦ "tiA.J] - »CjA.l]>/2. "cl.out"! 


/* s implication */ 


m[$x,$y,l] :0 

m[$x, $y, $2] - m[$y,$x,$z] 


0 


/* independent of theta 1 */ 
/* symmetry V 


/* flexible mass matrix is not */ 

/* function of flexible coordinates */ 


Do[i,3,n,Do[j ,3,n,m[i, j,$zl :01 } 


Do[i, 2,Do[j , 2,Do[lc,3 , 2,sym] | ] 

Do[i, 2 ,DoCj, 2 ,DoCk, 3 ,n,sym]] 

Do[i, 2 ,DoCj ,3 ,n,Do[lc ,3 ,n,sym] ] j 
Do[i, 3 ,n,DoCj , 2 ,Do[k, j , 2 ,sym] ] ] 

Do[i, 3 ,n, Do[ j , 2 , Do[fc, 3 ,n, sym] ] ] 


Figure D.4: Symbolic program of symmetric Christoffel symbol 


/. Rmflmceivm coupling fropmrty ot christoKml Symbol*/ 
l b-i, Put[c[j,i,W-=0 >>=■ i] •' CJ - oat * 1 1 

Do[i,2,Do[j,2,Do[]C,j, 2, refill 

DoCi, 2 ,DoCj, 2 , DoC3c,3,n, ref]]] 
DoCi,2 f DoCj,3 f nr D °CJ»J»n^jn] 

Do[i,3,n,Do[j ,2, Dock, j ,2, refill 
Do[i,3 ,n,Do[j ,2,Do[lc,3,n,refll 1 


Figure D.5: Symbolic program 


of reflective Christoffel symbol 
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/ 


* Output of Mass matrices Symbolic Program * 


'/ 


C2 

s2 

ul 

ule 

u2 


Cos[tii2] 

S in [ th.2 ] 

pll qll + pl2 ql2 
pile qll + pl2e ql2 
p21 q21 + ?22 q22 


/* Mass matrices */ 


mil • 


ml2 : 

m!3 : 
ml 4 : 
mlS : 
mlS : 
m2 2 : 
m2 3 : 
m2 4 : 
m2 5 : 
m2 6 : 
m2 3 : 
m3 4 : 
m3 5 : 
m3 6 : 
m44 : 
m45 : 
m4 6 : 
mSS : 
mS6 : 
m66 ; 


: :&3SSS3i 

al 4 l 2 -*l - a 2 *ll*pU« * 

a2*p21*x2 + a 2 *c 2 *ll*p 21 + a 2 ^p 2 is 2 

a2-p22-a2 * a 2 *c 2 *ll*p 22 * a2*p22-s2*ula 

a 2 *u 2 * 2 +- a2'*'X2 / '2 
a2*c2*plle*x2 + -(a2*pll«*s2 *u 2 ) 
a2*c2*pl2e*Tc2 + -(a2*pl2e*s2*u2) 
a2*p21*x2 
a2*p22*x2 

al*pll *2 + a 2 *plla ''2 

al*pll*p!2 +■ a2*plle*Pl2a 

a2*c2*plle*p21 

a2*c2*plle*p22 

al*pl 2^2 + a2*pl2e*2 

a2*c2*pl2e*p21 

a 2 * c 2 *p 12 e*p 2 2 

a 2 *p 21~2 

a 2 *p 21 *p 22 

a 2 *p 22 / '2 


m21:ml2 

m31:ml3; m32:m23 

m4l:al4; m42:m24; m43:m34 

m51:ml5 ; m52:m2S; m53:m35; aS4:a45 

m61:mlS ; m62:m2S; m63:m36; m64:m46, m65.m5 


Figure D.6: Symbolic program of velocity coupling matrices 
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/* Elements of Velocity Coupling Matrices V 
/•* Using Simplified Christoff el Symbol / 


c!21 

C221 
C131 
Cl4 1 

clSl 

C161 

C231 

C241 

C251 

C261 

C341 

C351 

C361 

C451 

C461 

c561 


ExCl/2*0(mll,th2] ] 

D[ml2,th2] 

Ex(l/2*D[mll,qll] ] 

Ex[ l/2*D[mll, ql2 ] ] 

Ex[l/2*D[mll,q21] ] 

Ex[l/2*D[mll,q22] ] 

1/2 (D[ml2 , qll] D[ml3,th2]) 

1/2 (D[ml2 ,q!2] D[ml4,th2]) 

Ex [ 1/ 2 ( D [ ml2 , q2 1 ] + DfmlS, th2] 
Ex [ 1/ 2 ( D [ ml2 , q22 ] +• D(ml6,th2]) 
1/2 (D[ml3 ,ql2] * D[ml4,qll]) 

1/2 (D[ml3 , q21] + D[ml5,qll]) 

1/2 (D[ml3 , q22 j + D[ml6,qll]) 
l/2(D(ml4,q21] *■ DCmlS, ql2]) 

1/2 (Df ml4 , q22 j +• D[ml6 , ql2] } 

1/2 (D[ml5 , q22 ] f D[ml6,q21]) 


] 

] 


C222 

Cl32 

C142 

C152 

C162 

C232 

C242 

C252 

C262 

C342 

C352 

C362 

C4S2 

C462 

C562 


1 / 2 *D [ m2 2 , th2 ] 
Ex[l/2(-D[ml3,th2 
Ex[l/2 (-D[ml4,th2 
l/2(-D[ml5,th2] f 
1/2 (-D[ml6 , th2 ] 
l/2*D(m22 , qll] 
l/2*D[m22,ql2] 
l/2*D[m22,q21] 
l/2*D[m22,q22] 

1/2 (D[m23 , ql2 ] +• 
1/2 (D[m23 ,q21] +• 
1/2 (D(m23 , q22] + 
l/2(D[m24,q21] 
!/2(D[m24,q22] +■ 
1/2 (D[m25 ,q22] + 


] + 0[m21, qll] ) ] 
] + D[m21,ql2])] 
D[m21,q21] ) 
D(m21,q22] ) 


D[m24 , qll] } 
D[m25 , qll] ) 
D[m26, qll] ) 
D[m25 ,ql2] ) 
D[m26,ql2] ) 
nr-rOli . <1211 1 


C223 

cl43 

C153 

C163 

C243 

C253 

C263 

C224 

Cl54 

Cl64 

C254 

C264 


Ex[l/2(-D[m22,qll] 
1/2 (-D(ml4 , qll] +■ 
1/2 (-D(ml5 , qll] * 
1/2 (-0[ml6,qll] + 
1/2 (-D(m24 ,qll] +• 
1/2 (-D[m25 , qll] +> 
1/2 (-D[m26 , qll] * 


+■ 2D[a32 , th2] ) ] 
D[m31,ql2] ) 
D[a31,q21]) 
D[m31,q22] ) 
0[m32,ql2] ) 

D(m32 ,q21] ) 
D[m32,q22]) 


Ex[l/2(-D(m22,ql2] 
l/2(-D(ml5,ql2] +> 
l/ 2 (- 0 [ml 6 ,ql 2 ] * 
1 / 2 (-D[m2S , ql2] +• 
l/2(-D[m26,ql2] +• 


+ 2D[m42,th2]) ] 
D[m41,q21] ) 
D(m41,q22] ) 
D[m42,q21]) 
D[m42,q22] ) 


C165 

C265 

C225 


1/2 (-D(ml6,q21] 
l/ 2 (- 0 [» 26 ,q 21 ] 
l/2(-D[a22,q21] 


+ D[a51,q22]) 
D(m52,q22] ) 

+• 2D[m52,th2]) 


C226 : 1/2 (-D[a22,q22] +■ 2D(m62,th2]) 


IS 

TV 


Figure D.6: (Continued) 
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appendix e 

Modeling of RALF using TREETOPS 

In this appendix, the modeling procedures of RALF using TREETOPS are 
described. 

TREETOPS requires several model definition data as follows: 

• Bodies 

• Modal Data 

• Hinges 

• Sensors 

• Actuators 

• Controllers 

• Function Generators 

• Interconnects 

• Devices (Spring/Damper) 

• Gravity 


• Constraints 
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The USER’S MANUAL for TREETOPS [73] is recommended reading for more 
detailed information. 

RALF can be expressed by six bodies which axe connected by six hinges as 
shown in Fig. 5.1. Bodies 1, 3. and 5 are assumed as rigid while bodies 2, 4. 
and 6 axe assumed as flexible. The modal data which must be supplied for each 
of the individual flexible bodies are denned in Table 5.1. These modal data can 

Table E.l: Modal data 


Modal Mass 

F. V'p 

Modal Stiffness 

F/OP 

Modal Damping (Optional) 

Fl?* 

Mode Shape of End Point 

0(0 

Mode Slope of End Point 

0'(O 

Modal Linear Momentum 

/ ipdm 

Modal Angular Momentum 

/ xtpdm 


be obtained from a NASTRAN output file. However, a NASTRAN/TREETOPS 
interface program is not available now. The MSC/PAL program can also provide 
the modal data. The mass properties and center of mass (C.o.M) of rigid bodies 
and the modal data of the first and the second mode of flexible bodies are shown 
in Table 5.2. TREETOPS uses MKS units. 

Note that the nodal mass matrix, N, is the unaugmented nodal masses and 
the mode shape matrix, 'P, is based on the augmented mode shapes. From the 
loaded interface component mode synthesis, only the augmented masses axe ob- 
tained. Therefore, the modal mass cannot be obtained from MSC/PAL directly. 
The modal linear momentum and the modal angular momentum also cannot be 
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Table E.2: Body input data and modal data of RALF 


Body Number 

Type 

Mass[Ivg] 

Inertia [N-m-s 2 ] 
C.o.M [m] 

[N/m] 

2 [N/m] 

E[=i V'.i m t /m T 

mi/ m T 

E[=i xtbiimi [Kg-m] 
E,=i xV\ 2 m i [Kg-m] 

0n(O 

i/»i2 (0 

*u(0 


1 

2 

Rigid 

Flexible 

9.9577 

8.3534 

2.3323 

13.5980 

0.4191 

1.1049 


7.224 


0.762 


2.2136E4 


43.0465E4 


0.2365 


-0.2886 


2.9148 


2.6581 


0.6054 


-0.0459 


-0.4043 


0.7S54 


3 

4 

Rigid 

Flexible 

5.9544 

8.3470 

0.7377 

18.944 

0.3048 

1.3208 


0.8052 


0.7S79 


0.7793E4 


2S.9S37E4 


0.2640 


-0.1408 


3.SS67 


-0.5129 


0.62S2 


0.6154 


-0.3262 


-1.1028 


5 6 

Rigid Flexible 

0.963 5.4465 

0.1193 S.S665 

0.3048 1.1049 

0.4804 
0.6187 
1.6711E4 
34.4127E4 
0.300 
0 

1.7151 

0.7611 

0 

0 

-0.5579 

-1.0134 


*«(0 
Notes: 

1) ^ r l , ; the node displacement vectors of 1st and 2nd mode. 

2) t^iit 0«'2 : the node displacement of node i of 1st and 2nd mode. 



199 


calculated by MSC/PAL directly. Therefore, a post processing program is required. 
The modal Unear momentum and the modal angular momentum are computed from 
the lumped mass method. In TREETOPS data, the modal linear momentum is di- 
vided by the total mass m T ■ The modal angular momentum terms are obtained 
from the cross product of the beam length axis and the Unk deflection. Therefore, 
the direction of the modal angular momentum is perpendicular to the plane defined 
by the x axis and the Unk deflection. Tachometer, position sensors, and torque mo- 
tors are located at the hinge 1 and 5. The decentralized PD controllers which are 
composed of the feedback of angular positions and angular velocities are added at 
the hinge 1 and 5. Interconnection between sensors, function generator, controller, 
and actuator are shown in Fig. E.2. The constraint is imposed between the body 3 
and 6. 

There are two basic steps to running the program. TREESET is an interactive 
preprocessor to define model data. It generates the PROBLEM.INT file. TREE- 
TOPS is a batch processor to run the PROBLEM.INT file. Time response is saved 
in the PROBLEM.PLT file. Mass, stiffiiess, and damping matrices are saved in the 
PROBLEM. AUX file. The Unear coefficients in matrix form (A,B,C,D) are saved 
in PROBLEM.MAT. 



VEL 1 
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Fig. E.2: Interconnection between sensors, controllers, and actuators 













201 


Bibliography 


[1] Asada, H., and Youcef-Toumi, K., “Analysis and Design of a Direct-Drive Arm 
with a Five-Bar-Link Parallel Drive Mechanism.” ASME Journal of Dynamic 
Systems, Measurement, and Control, 106, Sep. 1984, pp 225-230. 

[2] Asada, H., and Slotine. J.-J.E., “Robot Analysis and Control,” John Wiley & 
Sons, Inc , 1986. 

[3] Barbieri, E., and Ozgilner,U ., ” Unconst rained and Constrained Mode Ex- 
pansions for a Felxible Slewing Link,” ASME Journal of Dynamic Systems, 
Measurement, and Control, Vol. 110, Dec. 1988, pp 416-421. 

[4] Baumgarte, J., “Stabilization of Constraints and Integrals of Motion, Com- 
putational Methods m Applied Mechanics in Engineering, Vol.l, 1982, pp 1-16. 

[5] Bayo, E., “Computed Torque For The Position Control Of Open-Chain Flexible 
Robots,” Proc. IEEE International Conf. on Robotics and Automation, 1988, 
pp 316-321. 

[6] Benfield, W.A., and Hruda, R.F., “Vibration Analysis of Structures by Com- 
ponent Mode Substitution,” AIAA Jounal, 9(7), Jui. 1971, pp 1255-1261. 

[7] Book, W.J., Maizza-Neto, 0., and Whitney, D.E., “Feedback Control of Two 
Beam, Two Joint Systems With Distributed Flexibility,” ASME Journal of 
Dynamic Systems, Measurement, and Control, 97(4), Dec. 1975, pp 424-431. 

[8] Book, - W.J., “ Analysis of Massless Elastic Chains With Servo Controlled 
Joints,” ASME Journal of Dynamic Systems, Measurement, and Control, 
101(3), Sep. 1979, pp 187-192. 

[9] Book, W.J., Majette, M., “Controller Design for Flexible, Distributed Parame- 
ter Mech ani cal Arms Via Combined State Space and Frequency Domain Tech- 
nique,” Journal of Dynamic Systems, Measurement, and Control, 105, Dec. 
1983, pp 245-254. 

[10] Book, W.J., “Recursive Lagrangian Dynamics of Flexible Manipulator Arms,” 
The International Journal of Robotics Research. 3(3), Fall 1984, pp 8 <-101. 



202 


[11] Book, W.J., “Modeling, Design, and Control of Flexible Manipulator Arms: 
Status and Trends,” Proc. the NASA conference on Space Telerobotics, Jan. 31 
- Feb. 2, 1989, Pasadena, CA. 

[12] Boyle, J-M. Ford, M.P., and Maciejowski. J.M., “A Multivariable Toolbox for 
use with Matlab,” American Control Conference , June 1988, Atlanta. GA, pp 
707-712. 

[13] Burdick, J., “An Algorithm for Generation of Efficient Manipulator Dynamic 
Equations.” Proc. IEEE International Conf. on Robotics and Automation, San 
Francisco, CA, 1986, pp 212-218. 

[14] Cetinkunt, S., “Symbolic Modelling and Dynamic Analysis of Flexible Manipu- 
lators,” Proc. IEEE International Conf. on Robotics and Automation, Rayleigh, 
NC, 1987, pp 2074 - 2080. 

[15] Ceti nk unt, S.. “On Motion Planning and Control cf Multi-Link Lightweight 
Robotic Manipulators,” Ph.D. Thesis, School of Mechanical Engineering, 
Georgia Institute of Technology , Dec. 1987. 

[16] Chang, C.O., and Nikravesh, P.E., “An Adaptive Constraint Violation Stabi- 
lization Method for Dynamic Analysis of Mechanical Systems", ASME Journal 
of Mechanisms, Transmissions, and Automation in Design, 107, Dec 1985, pp 
488-493. 

[17] Chodas, J.L, and Man, G.K., “Design of the Galieo Scan Platform Control,” 
Journal of Guidance, Control, and Dynamics, Vol. 7, July - Aug. 1984, p422. 

[18] Chung, Y.C., Book,W.J., and Lu,S.S., “Modeling and Optimal Control of 
Lightweight Bracing Manipulator,” Internal Report, School of Mechanical En- 
gineering, Georgia Tech, May, 1986. 

[19] Craig, R.R. Jr., and Bampton, M.C.C., “Coupling of Substructures for Dy- 
namic Analysis,” AIAA Journal, 6(7), Jul. 1968, pp 1313-1319. 

[20] Doyle, J.C., and Stein, G., “Multivariable Feedback Design: Concepts for a 
Classical/Modem Synthesis,” IEEE Transaction on Automatic Control, 26, 
Feb. 1981, pp 4-16. * 

[21] Golla, D.F., Garg, S.C., and Hughes, P.C., “Linear State-Feedback Control of 
Manipulators,” Mechanism & Machine Theory, Vol. 16, 1981, pp 93 - 103. 

[22] Greewood, D.T., “Classical Mechanics”, Prentice - Sail, 1977. 

[23] Hablani, H.B.. “Constrained and Unconstrained Models: Some Modeling As- 
pect of Flexible Spacecraft,” Journal of Guidance and Control, 5(2), 1982, pp 
164-173. 



203 


[24] Hanafusa. H., Asada, H., and Mikoshi, T., “Design of Electrohydraulic Servo 
Systems for Articulated Robot Arm Control;’ IF AC Pneumatic & Hyrauhc 
Components , Warsaw. Poland, 1980. pp 223-228. 

[25] Hastings. G.G., “Controlling Flexible Manipulators, An Experimental Investi- 
gation ” Ph.D. Thesis. School of Mechanical Engineering, Georgia Institute of 

O 1 

Technology , Aug. 1986. 

[26] Hollerbach, J.M.. “A Recursive Lagrangian Formulation of Manipulator Dy- 
namics and a Comparative Study of Dynamic Formulation Complexity,” IEEE 
Trans. Systems. Man. and Cybernetics , 10(11), Nov. 1980, pp 730- (36. 

[27] Hou, S.N., “Review of Modal Synthesis Techniques and a New Approach, 
Shock and Vibration Bulletin , 40(4), 1969, pp 25-30. 

[28] Huggins. J.D., “ Experimental Verification of a Model of a Two-Link Flexi- 
ble, Lightweight Manipulator,” M.S. Thesis. School of Mechanical Engineering, 
Georgia Institute of Technology , June 1988. 

[29] Hughes, P.C., “Modal Identities for Elastic Bodies with Application to Vehicle 
Dynamics and Control,” Journal of Applied Mechanics . 47(1), Mar. I960, pp 
177-184. 

[30] Hughes, P.C., “Space Structure Vibration Modes: How Many Exist? Which 
Ones Are Important?,” IEEE Control Systems Magazine. FEb. 198 1 , pp 22-28. 

[31] Hurty, W.C., “Vibrations of Structural Systems by Component Mode Synthe- 
sis,” Proc. of the American Society of Civil Engineering , Aug. 1960 Vol86, pp 
51-69. 

[32] Ider, S.K., Amirouche, F.M.L., “Coordinate Reduction in the Dynamics of 
Constrained Multibody Systems - A New Approach, ’ Journal of Applied Me- 
chanics, Dec. 1988, Vol. 55, pp 899-904. 

[33] IMSL .Reference Manual, IMSL Inc, 1985. 

[34] Kane, T.R., and Levinson, D.A., “The Use of Kane’s Dynamical Equations in 
Robotics,” The International Journal of Robotics Research, 2(3), Fall 1983, pp 

3-20. 

[35] Kanoh, H., and Lee, H.G., “Vibration Control of One-Link Flexible Arm, 
Proc. of 24 th IEEE Conf. on Decision and Control, 3, Dec. 1985, pp H<2- 
1177. 

[36] Katie, D. and Vukobratovic, M., “The Influence of Actuator Model Complexity 
on Control Synthesis for High Performance RobotJTrajectory Tracking, I FAC 
Theory of Robotics , Vienna, Austria. 1986, pp 21 (-222. 


204 


[37] Kim, S.S., Vanderploeg, M.J., “QR Decomposition for State Space Represen- 
tation of Constrained Mechamcal Dynamic Systems. ’ ASME Journal of Mech- 
anisms, Transmissions, and Automation in Design , 108, Jun 1986, pp 183-188 

[38] Ivlema, V.C., and Laub, A.J., “The Singular Value Decomposition; Its Com- 
putation and Some Applications,” IEEE Transaction on Automatic Control , 
25, 1980, pp 164-176. 

[39] Liang, C.G., and Lance. G.M., “A Differential Null Space Method for Con- 
strained Dynamic Analysis”, ASME Journal of Mechanisms, Transmission, 
and Automation in Design , 109, Dec 1987, pp 405-411. 

[40] Krishnan, H., and Vidyasagar, M.. “Bounded Input Discrete-Time Control of A 
Flexible Be am : Theory, Simulation and Experimental Results , ASME Winter 
Meeting, Chicago, 1988, pp 1 - 16. 

[41] Low, K.H.. “A Systematic Formulation of Dynamic Equations for Robot Ma- 
nipulators with Elastic link,” Journal of Robotic Systems, 4(3) 198 i, pp 435- 
456. 

[42] Luh, J.Y.S., Walker, M.W., and Paul, R.P.C., “ On-Line Computational 
Scheme for Mechanical Manipulators,’’ ASME Journal of Dynamic Systems, 
Measurement, and Control, 102, June 1980, pp 69-76. 

[43] Luh, J.Y.S., and Zheng, Y.F., “Computation of Input Generalized Forces for 
Robots with Closed Kinematic Chain Mechanism, ’ IEEE Journal of Robotics 
and Automation , 1(2), Jun. 1985, pp 95-103. 

[44] MacFarlane, A.G.J. ed., “Frequency-Response Methods in Control Systems,” 
IEEE Press, 1979. 

[45] Maizza - Neto, 0., “Modal Analysis and Control of Flexible Manipulator 
Arms,” PhD Thesis, Dept, of Mechanical Engineering, MIT, 1974. 

[46] Mani, N.K., Haug, E.J, and Atkinson, K.E., “Application of Singular Value 
Deconposition for Analysis of Mechanical System Dynamics,” ASME Journal 
of Mechanisms, Transmissions, and Automation in Design, 107, Mar. 1985, 
pp. 82-87. 

[47] Meldrum, D.R., and Balas, M.J., “Application of Model Reference Adaptive 
Control to a Flexible Remote Manipulator Arm,” Proc. American Control Con- 
ference, 1986, pp 825-832. 

[48] Megahed, S., and Renaud, M., “Dynamic Modelling of Robot Manipulators 
Containing Closed Kinematic Chains,” Advanced Software in Robotics , Liege, 
Belgium, May 1984, pp 147-159. 

[49] Merritt, H.E., “Hydraulic Control Sytem,” John Wiley & Sons, Inc, 1967. 



205 


[50] Moon, F.C., “ Chaotic Vibration: An Introduction for Applied Scientists and 
Engineers,” John Wiely & Sons , 1987. 

[51] MSC/PAL2 Users Manual, MSC Inc, 1987. 

[5*5] Murray, J.J., and Neuman, .C.P.. ‘‘ARM: An Algebraic Robot Dynamic Mod- 
eling Program.” Proc. of the First International IEEE Conference m Robotics, 
Mar. 1984, pp 103-114. 

[53] Nelson, W.L., and Mitra, D., “Load Estimation and Load Adaptive Optimal 
Control for a Flexible Robot Arm,” Proc. IEEE International Conf. on Robotics 
and Automation , San Francisco, CA, Mar. 1986, pp 206-211. 

[541 Nikravesh, P.E, “Some Methods for Dynamic Analysis of Constrained Mechan- 
ical System: A Survey,” Computer Aided Analysis and Optimization of Me- 
chamcal System Dynamics, E.J.Saug, ed., Springer- Verlag, Heidelberg, West 
Germany, 1984. 

[55] Ogata, K, “Modem Control Engineering,” Prentice- Hall, 1970. 

[56] Ower, J.C., and Van De Vegte, J., “Classical Control Design and a Flexible 
Manipulator: Modeling and Control System Design”, IEEE Journal of Robotics 
and Automation , RA-3, 5, Oct., 198 <, pp 485*489. 

[571 Park, T.W., and Haug. E.J., “A Hybrid Numerical Integration Method for 
M achin e Dynamic Simulation”, ASME Journal of Mechanisms, Transmissions, 
and Automation in Design , 108, Dec 1986, pp 211-216. 

[58] Pfeiffer, F., Gebler, B., and Kleemann, U., “On Dynamics and Control of 
Elastic Robots,” IFA C -Symposium on Robot Control , 1988, pp 4. 1-4.5. 

[591 Punyapas, N., “The design of Hydraulic Valve Circuit under^ Microprocessor 
Control for Increasing efficiency of Hydraulic Servo System, Ph.D. Thesis, 
School of Mechanical Engineering, Georgia Institute of Technology , Mar. 1982. 

[60] Rovner, D.M., and Cannon, R.H.,Jr., “Experiments Toward On-line Identifica- 
tion and Control of a Very Flexible One-Link Manipulator, The International 
Journal of Robotic Research, 6(4), 1987, pp 3-19. 

[611 Sakawa, Y„ Mauno, F., and Fukushima, S., “Modeling and Feedback Control 
of a Flexible Arm,” Journal of Robotic Systems, 2(4), 1985, pp 453-47- 

[62] Schmitz, E., “Experiments on The End-Point Position Control of a Very Flex- 
ible One-Link Manipulator,” Ph.D. Thesis, Department of Aeronautics and 
Astronautics, Stanford University, Jun. 1985. 

[63] Schutter, J.D., et al., “Control of Flexible Robots Using Generalized Nonlinear 
Decoupling,” IFAC-Symposium on Robot Control, 1988, pp 98.1-98.6. 


206 


[64] Siciliano, B. and Book, W.J., “A Singular Perturbation Approach to Control 
of Lightweight Flexible Manipulator,” The international J ournal of Robotics 
Research . 7(4), 1988, pp 79-89. 

[65] Siciliano, B., Yuan, B., Book, W.J., “Model Reference Adaptive Control of a 
One Link Flexible Arm,” IEEE Conf. Decision and Control , 1986, pp 91-95. 

[66] Shabana, A. A., “Substructure Synthesis Method for Dynamic Analysis of 
Multi-Body Systems,” Journal of computers and Structures , 20(4), 1985, pp 
737-744. 

[67] Singh, R.P., and Linkins, P.W., “Singular Value Decomposition for Constrained 
Dynamic Systems,” ASME Journal of Applied Mechanics , 52(4), Dec. 1985, pp 
943- 948. 

[68] Skelton, R.E.. Hughes, P.C., and Hablani, H.B., “Order Reduction for Models 
of Space Structures Using Modal Cost Analysis,” Journal of Guidance and 
Control , 5(4), 1982, pp 351-357. 

[69] SMP Reference Manual. Inference Corporation , 1983. 

[70] Sunada, W., and Dubowsky, S., “On the Dynamic Analysis and Behavior of In- 
dustrial Robotic Manipulator with Elastic Members,” Journal of Mechanisms, 
Transmissions, and Automation in Design , 105, Mar. 1983, pp 42-51. 

[71] Tourassis, V.D., and Ne uman , C.P., “Properties and Structure of Dynamic 
Robot Models for Control Engineering Applications,” Mechanism and Machine 
Theory , 20(1), 1985, pp 27-40. 

[72] Tourassis, V.D., and Neuman, C.P., “The Inertial Characteristics of Dynamic 
Robot Models,” Mechanism and Machine Theory , 20(1), 1985, pp 41-52. 

[73] TREETOPS Users Manual, DYNACS, Clearwater, Florida, 1989. 

[74] Tsujisawa, T., “A Reduced Order Model Derivation for Lightweight Arms with 
A Parallel Mechanism,” Proc. IEEE International Conf. on Robotics and Au- 
tomation, Scottsdale, AZ, May 1989. 

[75] Venkatr aman , V. and Mayne, R.W., “Dynamic Behavior of a Hydraulically 
Actuated Mechanism part 1: Small Perturbations,” Journal of Mechanisms, 
Transmissions, and Automation in Design , 108, June 1986, pp 245-249. 

[76] Wehage, R.A., and Haug, E.J., “Generalized Coordinate Partitioning for Di- 
mension Reduction in Analysis of Constrained Dynamic Systems,” ASME 
Journal of Mechanical Design , 104(1), 1982, pp 247-255. 

[77] Wie, B., Lehner, J.A., and Plescia, C.T., “Roll/Yaw Control of a Flexible 
Spacecraft Using Skewed Bias Momentum Wheels,” Journal of Guidance, Con- 
trol, and Dynamics, Vol. 8, July-Aug. 1985, pp 447-453. 



207 


[78] Wie, B., and Byun, K., “New Generalized Structural Filtering Concept for 
Active Vibration Control Synthesis,” Journal of Guidance and Control, 12 (2), 
Max.-Apr. 1989, pp 147-154. 

[79] Wilson, T.R., “The Design and Construction of a Flexible Manipulator,” M.S. 
Thesis, School of Mechanical Engineering, Georgia Institute of Technology , 
Dec. 1985. 

[80] Yang, Y.P., and Gibson, J.S., “Adaptive Control of a Manipulator with a 
Flexible Link,” Journal of Robotic Systems, 6(3), 1989, pp 217-232. 

[81] Yocum, J.F., and Slater, L.I., “Control System Design in the Presence of Se- 
vere Structural Dynamics Interactions,” Journal of Guidance, Control, and 
Dynamics, Vol. 1, Mar.- Apr., 1978, pp 109-116. 

[82] Yuan, B. -S., “Adaptive Strategies for Position and Force Controls of Flexible 
Ar ms ,” Ph.D. Thesis, School of Mechanical Engineering, Georgia Institute of 
Technology, Mar. 1989. 

[83] Yuh. J., “Application of Discrete-Time Model Reference Adaptive Control to 
a Single- Link Flexible Robot,” Journal of Robotic Systems, 4(5), 1987, pp 621- 
630. 

[84] Yurkovich, S., and Pacheco, F.E., “On Controller Timing for a Flexible-Link 
Manipulator with Varying Payload,” Journal of Robotic Systems, 6(3), 1989, 
pp 233-254. 



208 


VITA 

Jeh Won Lee was bom in on He received the 5.S. 

degree in Mechanical Engineering from Seoul National University in 1979 and the 
M.S. degree in Mechanical Engineering from Korea Advanced Institute of Science 
and Technology (KAIST) in 1981. After graduation, he worked in Daewoo Motor 
Company for one year and in Daewoo Heavy Industries for three and half years. 

During that time, he was married to a lovely lady, Kyung^^l and had one 
daughter, Yoo Jin. 

He started his doctor degree program at the school of mechanical engineering of 
Georgia Institute of Technology in January, 1985. He was a teaching assistant and 
a research assistant. 

His master thesis topic was An application of PUD control law to machine tool 
chatter control using electro - hydraulic actuator. In Daewoo Heavy Industries, he 
was a key member of development of an arc welding robot. He currently works 
at Pointing Control Branch of Structures & Dynamics Division of Marshall Space 
Flight Center in Huntsville. His current interesting research area is Robotics, Dy- 
namics and Control of multi-body flexible structures. 




